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Similarity Laws for Stressing Heated Wings 


H. S. TSIEN* 
California Institute of Technology 


SUMMARY 


It wili be shown that the differential equations for a heated 
plate with large temperature gradient and for a similar plate at 
constant temperature can be made the same by a proper modi- 
fication of the thickness and the loading for the isothermal plate. 
This fact leads to the result that the stresses in the heated plate 
can be calculated from measured strains on the unheated plate by 
The applica- 


a series of relations, called the ‘‘similarity laws.” 


tion of this analog theory to solid wings under aerodynamic heat- 
ing is discussed in detail. The loading on the unheated analog 
wing is, however, complicated and involves the novel concept 
of feedback and ‘“‘body force’’ loading. The problem of stressing 
a heated box-wing structure can be solved by the same analog 


method and is briefly discussed. 


INTRODUCTION 


es HIGH STAGNATION TEMPERATURE for flight of 
aircraft at supersonic speeds results in severe aero- 
dynamic heating of the surfaces of aircraft. For in- 
stance, in a recent paper by Kaye,!' the transient tem- 
perature distribution in a wedge-shaped solid wing was 
calculated for accelerated flight and was found to be 
rapidly varying both with respect to the space points 
in the wing and with respect to time instants. When 
there are large temperature gradients in the material, 
there are generally large thermal stresses due to the 
uneven thermal expansion of the material. Therefore, 
concurrent with the severe aerodynamic heating, there 
is the problem of determining the thermal stresses in 
the wing. The purpose of this paper is to suggest a 
method of stressing such heated wings. 

The starting point of the stress analysis is the tem- 
perature distribution in the wing. Because the rate of 
change of temperature is small in comparison to the 
speed of sound propagation in the material, the stress 
calculation can be considered as a quasi-steady problem. 
That is, the stress in the wing at each time instant can 
be calculated from the temperature distribution at that 
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time instant without the inertia effects from the time 
variation of the material displacements required by the 
changing stresses and the thermal expansions. This 
problem in general theory of elasticity is greatly sim- 
plified by the fact that the wing is thin, and therefore 
the Kirchhoff hypothesis for bending of thin plates 
holds. In fact, the problem of thermal stresses in an 
elastic plate has been treated some time ago by Nadai.’ 
The present study generalizes the earlier theory in two 
directions: Nadai’s assumption of linear temperature 
profile across the plate is no longer necessary. The 
temperature profile is now arbitrary. Secondly, 
Young’s modulus EF of the material is now allowed to 
vary as a function of temperature. Thus, the effects of 
decrease in Young’s modulus with increase in temper- 
ature can be taken into account in the present 
theory. 

However, the main purpose of the present paper is 
not just to construct a theory for heated plates. The 
main purpose is to utilize the theory to formulate a law 
of similarity, whereby the problem of stressing a heated 
wing can be solved by performing a set of experiments 
on a properly proportioned and properly loaded wing 
at room temperature without heating. This concept 
of similarity law has been explored before by the pres- 
ent author*® for thin-walled cylinders. It is believed 
that this approach to the problem of stressing a heated 
wing has many practical advantages over a pure analyti- 
cal solution, as will be discussed later in the paper. 


Basic EQUATIONS FOR HEATED PLATE 


Let the plane of the plate be the x-y plane, and its 
variable thickness be b(x,y). The z-axis is normal to 
the plate and pointing downward. The temperature 
above the ambient is specified by T(x,y; 2). Actually, 
the distance z is measured from the ‘“‘median surface,” 
the position of which will be specified presently. There- 
fore, it is implied here that the median surface, al- 
though not exactly coinciding with the x-y plane, is 
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nevertheless sufficiently flat to be considered as the 
x-y plane. Let u, v, and w be the displacement of a 
point (x,y) on the median surface in the x, y, and z 
directions, respectively, due to elastic strain and thermal 
expansion. Then, according to Kirchhoff’s bending 
hypothesis, the total strains of a point (x, y, z) in the 
plate are given by 


er, = (Ou/dx) — 
2(0°w/ Oy") (1) 


€é, = (dv/dy) — 
(0u/Oy) + (Ov/Ox) — 22(0°w/Oxdy) 


2(0?w/Ox?) 


Tx = 


where e¢, is the direct strain in x, ¢€, is the direct strain 
in y, and y,, is the shear strain in the x-y plane. All 
other strains are small and negligible for thin plates. 

Let E(T) = E(x,y; 2) be the variable Young’s modu- 
lus, vy the constant Poisson’s ratio, and oz, oy, Tz, the 
significant components of stresses. Then the strains 
can be computed from these quantities as follows: 


;2)] + aT(x,9; 2) 
z)] + aT (x,y; 2) (2) 


ér = [(o, — vo,)/E(x,y 


¢, = [(e, — voz)/E(x,y; 


Tx = 2(1 + V) Try/E(x,y; z) 


where a is the coefficient of thermal expansion. By 
solving Eq. (2) for stresses, one has 
_E 3 yi 8 z) 


oc = 1 [(€z -+- vey) 


o 


— (1 + v)aT (x,y; 2 


E(x,y; 2) 3) 
se id [(ey + ver) — (1 + v)aT (x,y; 2 ( 


2)/2(1 + v)] Yry 


tx = [(E(x,9; 


In the theory of thin plates, the important quantities 
are not the stresses but are rather the sectional forces 
and the sectional moments derived from the stresses. 
These sectional quantities are defined as 


as S 0, dz, N, = Ps dz, Ny = S tr dz (4) 


= fo zdz, M, = fo,zdz, M, = 
—S try2 dz (5) 


where all integrations extend across the whole thickness 
of the plate. N, and N, are then sectional normal 
forces, N,, is the sectional shear, 1/7, and M, are the 
sectional bending moments, and M,, is the sectional 
twisting moment. By substituting Eq. (1) into Eq. 
(3) and then into Eqs. (4) and (5), one has, for example, 


V.=D (= 4 2") D (= 4 4) V 
sie = 0 ae Vv oy 1 ax? v dy? iV7 
Ou Ov O*w O*w 
M. = D —]—D, — Mr 
= ss =<) oe ial oA . 
where 
l a 
Do(x,y) = af E(x,y; 2) dz (6) 
= 
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I 
D(x,y) = = [Hess 2)2 dz (7 
l-— pv 
l , 
D.(x,y) = E(x,y; 2)2? dz (8 
l— vp ‘i 
; a on = 
Nr = E(x,y; 2)T (x,y; 2) dz (9) 
l—vyp 
a — so 
Mr ms E(x,y; 2)T (x,y; 2)2 dz (10 
> 


Ny and Mr being due to thermal expansion of the ma- 
terial can thus be called thermal sectional normal force 
and thermal bending moment, respectively. From the 
above typical expressions for NV, and M,, it is seen that 
a simplification results if the median surface is so chosen 
that 


D, = 0 (11 


In fact, this is the condition for fixing the position of 


the median surface. With the median surface so deter- 


mined, the sectional quantities can be calculated as 


Oo Ov : 
a im De ( - ) = | 


ox oy 
Ov Ou | 
N, =D -} — Nr ¢ 12 
, ? S bith Ox tf \ 
. gierte (or + ot) | 
oe ? Oy Oy 
O*w O*w ) 
M, = —D, — Mei 
es i =) . | 
0*w 0*w 
uw —p.( , *) —~Mr} (13 
oy" Ox" | 
o*w | 
M, = l— Dz | 
; ¥) Oxoy 


The equilibrium of forces in the medium surface re- 


quires 
(ON,/Ox) 
(ON,,/0x) 


+ (ON,,/dy) = 0) 


(14) 
+ (ON,/dy) = 0) 


By substituting Eq. (12) into Eq. (14), the following 
two equations for u and v are obtained: 


re) | (= x *) ‘ } 
Ox "XO: ” ov 
1-—v0 | (= ) ONr 
Do + = 
2 oy Ox oy Ox | 


is o | (= 4 on) 4 
2 CO ” Lox oy 
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SIMILARITY LAWS FOR STRESSING HEATED WINGS 3 
fied in terms of forces instead of displacements, it would 
be more convenient to use the stress function ¢(x,y) 
as the basis of solution. & is defined through the fol- 
lowing relations and satisfies Eq. (14) automatically: 


When the temperature distribution is specified, the 
right sides of these equations are known. Then Eq. 
(15) is a set of simultaneous partial differential equa- 
tions for u and v. If the boundary conditions of the 
problem are specified in the displacements, Eq. (15) . 
js the proper basis of solution. — 0*y a 0*y —— O° (16) 
If the boundary conditions of the problem are speci- Oy?’ ie Oxdy "Ox? 


By substituting these relations into Eq. (12), one has 


Q oy dy ] 
> . ( ——— ‘) ti Ne | | 
Ox (1 — v?)Do oy" Ox* | 


Ov 1 0° **) rte | 
dy (1 — v°)Do (& ” Oy? . atid ; 

Ov Ou 2 0*y 

i ( + ) = 2 

Ox Oy (1 — v)Do Oxdy J 


When the first equation of Eqs. (17) is differentiated twice with respect to y, the second equation is differentiated 
twice with respect to x, the third equation is differentiated with respect to x and y, and the results are added, a 


single equation for ¢g is obtained: 
1 oO? [ 1 /d’¢ "9 2 0” 1 Oy 
Ls eLACS--BN + tases) 
(1 — v?) Ox? LDo \Ox? oy? (1 — v) Oxdy \Do Oxdy 
l 0 1 /O°¢ 0’ l N 
o (oe? “)|+ v2 (7) = 0 (18) 
1 — pr Oy" Dy oy" Ox? ] + v Do 


where V’ is the Laplacian operator. At the boundary of the plate (Fig. 1a), if /, m are the direction cosines of the 
“outside’’ normal ” to the boundary direction s, then sectional normal force N,, to the boundary and the sectional 


shearing force V,,, along the boundary are given as follows: 


N,+N, N,—N, l | (= A) 20 
i, = - (]? — m? N,,(2lm) = — A’e (/? — m? - — 2lm —— 19) 
+ m*) + wf 5 ae + 5 m”) mer Dx? nap ( 


9 ») 


oO” “a oO” a oO” a 
Nis = —(N, — N,)lm + (? — m*)N,, = — lm ( s -- ‘) — (l? — m?) = (20) 
Oy? Ox? xoy 


With the boundary described and the boundary forces specified, Eqs. (18), (19), and (20) then give the proper 


formulation of the problem. 
If p is the downward pressure loading on the plate, then the equilibrium of forces normal to the plane of the plate 


requires 
0°?M, 0°M., . 7M, _ Ow . 0 , ow 
——2 — = —p — N,— — 2N,, - N= 
Ox’ OxOoy Oy" Ox* Oxoy Oy? 
By substituting the sectional moments obtained from Eq. (13) into the above equilibrium equation, an equation 


for the lateral deflection w results: 


0° |p: (= P- ~~) | Rial o? (o. = ) 4 0? | D: (= re 4) Z 
Ox? | ~ \dx? Oy? 7 Oxdy \ ~ Oxdy oyv?L ~~ \Ody? Ox 

4.9N O*w : o*w 

ane Saedy Dy? 


O*w 
—-V?M,y+p4+N, ax? (21) 


x 
Part of the boundary conditions will be expressed in terms of w directly. For instance, if the boundary s with 
outside normal 7 is a “‘fixed edge,’’ then on that boundary 


(22) 


w= 0, ow/on = 0 
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The other part of the boundary conditions is described in terms of sectional moments. 


cosines of the outside normal, the sectional bending mc 
the boundary (Fig. 1b) are, then, 


ICAL SCIENCES—JANUARY, 1953 


If /, m are the direction 
ymments M,, M, and the sectional twisting moment J/,, on 


M, M M, — M l 
M, = - ke = + ——* (I? — m*) — Myz,(2lm) = —My — " DV — 
1 — Ow dow *ay) 
. ” DL? — m’) ( ie *) — 2(1 — v)Delm 
2 Ox? Oy? Oxoy 
t (23) 
VM, M M, — M, l 
M, us . (7? — m*) + M,, (2lm) = —Mp — . 4 DNV?*w + 
l—yp Ow dw 072 
D,(l? — m?) _ 2(1 — v)Dz 
, ( m (= -) + 2( v)Dolm — 
" , Ow dw 2 ,. Ow 
M,, = (M, — M,)lm + (P -— m’*)M, = -(1 - Dams - dy? ‘) + (1 — v)D,(l? — m?) dxdy (24) 


For a “free edge,’’ the Kirchhoff boundary conditions 


are, then, 
OM, M, — M, OMas 
Lh Sn 


= 95 
On r Os o @) 


M, = 0, 


where 7 is the radius of curvature of the boundary s, or 


1/r = (1/m) (dl/ds) (26) 


PLATE AT CONSTANT TEMPERATURE 


For a plate at the constant reference room temper- 
ature, all the temperature terms in the equations de- 








(2) 
Ac > x 
Pst 
Nas 
cn nN, 
y Z 
(6) 










ame ges 


°™, 


Fic. 1. 


Boundary forces and moments. 


veloped in the previous section vanish. By using a 
bar over a quantity to differentiate the present case 


from the previous case, one has for strains 


é, = (0%/d0x) — 3(0°/dx") ) 


éy = (00/O0y) — 2(0°@/dy?) (27) 
Vry = (O00/0x) + (O%/Oy) — 23(0%/Oxdy) J 
Let E be Young’s modulus, now a constant, then 
a; = [E/(1 — v)](&+» vey) | 
oy = [B/(1 — v*)] Cy + vez)? (28) 
Try = [E 2(1 + v) Vey ) 


The sectional forces and the sectional moments defined 
as Eqs. (4) and (5) are 


a a 4 o 
N,; = Do 








N, = Dy 29) 
| e* , ” 
a 1 edhe D es 22) 
- "Nox * dy/ j 
a, = —D 3 nd 22) 
Nant "Oy? | 
i D ( aN -) | 
See 2 Begg f 30) 
' Oy? Gn} | 
. «1 «sh | 
—s Me? axdy } 


where, with b(x,y) as the thickness of the plate at the 
point (x,y), 


- b(x,y), D.(x,y) = - 2 (31) 
—— 1—»? 12 


Do(x,y) = 
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The median surface now lies midway between the upper 0 7 Ou Od } 
and the lower surfaces of the plate. ox : ox a oy +r 
Now introduce a potential F(x,y) for the body ’ 
forces X and Y, which are the forces per unit 2a |D, (= =] = oF 
area Of the median surface in the x and y direc- 2 oy Ox oy Ox 
tions: . + (34) 
; l1—vd |p (5° 4 oY nq 
X = -dF/ox, Y= -—-dF/y (32) 7 ald 
ra) | D (5 ‘ *) OF 
Thea 2 ‘ ; . Fa 41; ; } rf eC > > ) Vv = 
Then the equations of equilibrium of forces in the me dy *\oy ae dy 


dium surface are 
5 m Eq. (33) will be automatically satisfied, if a new 
(ON,/Ox) + (ON7,/Oy) — (OF /ox) = O| (33) Stress function ¢(x,y) is introduced such that 
(ON,,,/0x) + (ON,/dy) — (OF/dy) = 0f 

' i - 0°¢ ~ 0°¢ - 0°¢ 


= -F, —Ny = » N= +F 
oy" . Oxoy Ox? r 


tel 
By substituting Eq. (29) into the above equations, a 
set of equations for # and 7 is obtained: (35) 


The equation for 2 is now 


I Oo? E (ee | nN 2 o? (; 2.) 4 
1— v2dx? LD) \dx? dy? 1 — vdxdy \Do dxdy 


1 oO | (o ~) , ( *) ; 
Sah pa v?(=)=0 (36) 
1 — v? Oy? LDo \Oy? Ox? ee ae Do 


At the boundary of the plate, 


- N, N, N, — N, re 
Kh = - oo = - (12 — m?) + N,,2lm 
1 o*¢ or) oe 
= 7 V’°o (J? — m?) — — 2lm (37) 
* 2 * 2 ee Ox" Oxoy ™ 
N,. = —(N, — N,)im + (P - m?) Nr 
O°o 0° 5 : 2D 
= —Ilm - — ——]} — (/? — m’*) (38) 
Oy" Ox? Oxoy 


If p(x,y) is the downward pressure loading of the plate at the point (x,y), the equilibrium of lateral forces gives 


the following equation for @: 


| D (= ' a4) 191 o? (p | 4 roe | D ( i 1 4] 
be . ) - 9 —— " x = 
Ox*L ~~ \ Ox? oy? ; Oxdy " Ody oy? ~ \oy? ” Ox? 
5+ 408,25 +h 
Ox? ‘ Oxdy ' Oy? 
The sectional bending moments M7,,, M7, and the sectional twisting moment M,,; are, then, 
= l = l—vp-;- Ow O"a = O*a 
fa ie ihe. of - ) — 2(1 — v)Delm 
2 2 Ox Oy" OxOoy 
(40) 
_ i ++ > l-—vy- "a O°a a 0°a 
M,=- ; DV *® D,(l? — m?*) -_— ) + 2(1 — v)Delm | 
2 . Ox? Oy" OxOoy | 
"a 0°w "a 
(41) 


- (1 — v)D,(l? — m*) 
Ox? =) r ii ” Oxoy 
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Loads and boundary conditions of a solid heated wing 


SIMILARITY LAWS FOR SOLID THIN WINGS 


By comparing the corresponding equations for the 
heated plate and for the plate at room temperature, 
their similarity is evident. The question is then 
whether it is possible to find a corresponding plate at 
room temperature such that, with proper loading, the 
plate will give solutions that are similar to that of the 


heated plate. To fix ideas, let the heated plate and the 


unheated plate both have the same plane form. The 

Oo? | ( =w fi =) +901 o” (p ow 
) + 2 ys = 2 

Ox? “ \Oe Oy" ; Oxdy ~ OxOY 

l 

Ap 


Therefore it can be deduced, by comparing the above 

equation with Eq. (21), that, in order for w and @ to 

satisfy the same equation, 
1=8 


p = —BNrV"® = BuV? My + Bup 


(48) 
(49) 


On the boundary of the heated plate, the forces JN, 


and N,, vanish. By comparing Eqs. (19) and (20) 
with Eqs. (37) and (88) and by noting Eqs. (43) and 


(44), 


N, ai BN7, Nas = 0 (50) 


Similarly, by comparing Eqs. (23) and (24) with Eqs. 
(37) and (388) and by noting Eqs. (46) and (47), the 
boundary conditions of Eq. (25) for the free edge of the 
heated plate reduces to 
MN, = BuM yr 
OM, M, — M, OM ns oMr\ 
> + = Bu ( 


(51) 
On r Os on 
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temperature distribution of the heated plate is speci- 
fied so that Nr and >» are known functions of the co. 
ordinates (x,y). The heated plate is built in at one end 
but free at all other edges (Fig. 2). The lateral pres. 
sure loading p for the heated plate is also specified, 
The problem is to find the thickness distribution b(x,y 
and the loading of the unheated plate, such that P(x,y 
and @(x,y) are the same as g(x,y) and w(x,y), respec- 
tively, except a proportionality factor. 

Let a and 8 be two constants; then it is clear that the 
Eqs. (18) and (36) for g(x,y) and ¢(x,y) will be exactly 
the same if 


Dy = aDy (42) 
% = Be (43 
F = BN; (44 
Then, because of Eqs. (16) and (35), 
N, = B(N, + Nr), Ny = BNy,, N, = 
B(N, + Nr) (45 
Now, if \ and uw are two additional constants, and 
Do. = XD» (46 
w= pw (47) 


then the equation for @, Eq. (39), can be rewritten as 
) 4 oO” | ( "Wy 4 | 
- Vv — 
Oy" Oy" ox 
b + B Vp 2@ + B (v “Ww Lon O°w LN “| 
y Nr V*@ N, + 2N, - N, 
mM d * ac? ' OXOY ' Oy? 


The boundary conditions for the built-in edge, Eqs. (22), 
give 


Ow/On = 0 (52) 


As a first application of the theory developed, con- 
sider the simple case of a thin wing of solid section: 
The purpose is to find the stresses in the heated wing 
by performing experiments on an unheated “‘correspond- 
ing’ or analog wing. The first step is to understand 
clearly the meaning of Eqs. (42) and (46). These 
equations, in fact, determine the thickness 4 of the 
Now let g be ratio E/E and 7 the non- 
the 
—1 at upper 
Let m be the 


unheated wing. 


dimensional thickness variable measured from 
upper surface of the plate such that » = 
1 at the lower surface. 


Then, 


kb “1 
a g(x,y; n) dn 
2( l a. ory i 1 


surface and 7 = 
value of n at the median surface. 
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i - With the geometry of the analog wing determined, 
Dixy) = 11 — y?) / g(x,¥5 0) (n — mo) dn the next step is to specify the loads on it. Eq. (50) 
ie shows that the unheated wing must be loaded with a 
a Eb? - sectional tension force BN7 at the boundary but no sec- 
ae tee D:(x,y) = ft = of / g(x,ys 0) (n — m0)* dn tional shearing force. The first of Eqs. (51) shows 
—e end that, at the free edge, the analog wing should have a 
ital By using Eq. (31) and the above relations, the condi- bending menuasee equal wea BuMr. ieee numerical 
b(x.y) tions of Eqs. (11), (42), and (46) can be written as value OI mw 1s arbitrary and a - the disposal ol the ex- 
i perimenter. The second of Eqs. (51) can be inter- 
?(x,y) “1 preted as an upward support force —8u(0.\/7/On) per 
es / g(x,y; 0) (n — mo) dn = 0 (53) unit length of boundary but no twisting moment at the 
wit boundary. Eq. (52) shows that a built-in edge of the 
stabs: © iw Os o oe . heated wing corresponds to a built-in edge of the analog 
<actly 5 O(%Y) / g(x,y; n) dn = b(x,y) (54) - 
. J -1 wing. 


\ - | On the unheated wing, Eqs. (32) and (44) specify a 
(42) ~ (x,y) Fj g(x,v; n) (n — mo)? dn = , b*(x,y) (55) body force 
S y : 


s eta laa : ee ee (57) 
By eliminating }, 6, and mo from these equations, one Ox Oy 
(44) te 
' — These are forces in the x and y direction per unit area 
| at - _ of the wing surface and are loads perhaps new to struc- 
(-, “) ( / g in) = (/ g-n° in ) x tures testing. The lateral pressure loading p on the 
z ‘ini es iia analog wing is specified by Eq. (49). If Az, and Az, 
ae “4 ‘ denote the differences of x and y components of direct 
( / gf in) a ( / g°9 in) (96) — strains measured on the isothermal wing on the top and 
= vale the bottom surfaces—i.e., according to Eq. (28) 
ae Eq. (56) shows that the ‘“‘Young’s modulus profiles’’ 5 5 : oa 
(47) g(n) for various points (x,y) of the heated plate are not Aé,(x,v) = iss on °) _ is y ) = b(x,y) : 


entirely arbitrary but must satisfy that relation so 
n as that the ‘“‘similarity’’ between the heated wing and the b b Ow 
isothermal analog wing is possible. If g is a fixed Aa,(x,y) = % (ss = ;) = Ss (2: 4) = b(x,y) y2 
function—+.e., the temperature profiles across the plate (aa 
are similar—or Young's modulus is constant, then Eq. 
(56) is certainly satisfied.* For any given problem of a — then Eq. (49) can be written as 
te pry wing a, aed first step is " compute the — b = (—BNr/B) (dé, + Az) — BuV?Mr + Bup (59) 
) of a*®/X\ or a®/B for various points of the wing by Eq. 
(56). If these values do not differ very much, then, as Needless to say, if there are concentrated lateral loads, 
an approximation, the average of the computed a*/8 the ratio of these loads on the test wing and on the origi- 
can be used. Then the similarity procedure is possible. nal wing is also equal to Bu. 
2) But in any event, the value of a*/@ is fixed by the prob- The physical significance of Eq. (59) is different from 
lem, not at the free choice of the stress analyst. Thus, that of Eq. (57) in that the body force loading is com- 
if Bis chosen, then a is fixed. Furthermore, \ is equal pletely specified before the test because Nr is a known 
52) to B according to Eq. (48). Therefore, out of the four quantity, while the first term of p itself depends upon 
constants, a, 8, A, and uw, only two (8 and uw) can be the experimentally determined Az, and Az,. There- 


chosen arbitrarily. When this is done, Eq. (54) deter- fore, if the loads are considered as input to the wing 


nl- A P ; = ‘i R ‘ 
i mines the appropriate thickness )(x,y) for the unheated structure and the experimentally measured strains as 

Tl. ; * _ pe ? 

Me analog wing. output, then Eq. (57) shows that the output also partly 
5 * . . . 

ae ’ determines the input. In other words, there is a 

d- * Mathematically speaking, Eq. (56) is a functional equation leter _ ert po ; : 

- The solution of it gives the explicit required character of g. This feedback link"’ in the experimental setup of the analog 
si. problem is discussed in the Appendix. wing. 

1e 

- 

le " é 

- With the loading on the unheated wing so specified, the relation between ¢ and ¢ and between @ and w are given 
e by Eqs. (43) and (47). The stresses in the original heated wing can then be computed from the measured strains 


on the test wing. Forinstance, from Eqs. (1) and (3), 


E(x,y; 2) [é 4 o) ( Ww ra =~) a+ T\ | 
(x,y; 3) = — 2 - \aT (x,y; 2 
oe 1 — vy? L\Ox , Oy Ox? , Oy’ ‘ 
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By using Eq. (17), the above equation is reduced to 
E(x,y; 2) {Nr 1 O° 


1—>* WD” De oy? 


o,(x,y; 2) = 
But according to Eqs. (42), (43), (29), and (35), 


1 0p a 1 0p a i 
Do Oy" 7 B Dy oy? 7 B Dy 


6 


x) 


However, 


Ou i? @a.. 
ox 2b“ ~ 9 ° 


[N, 
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O*w O*w = 
— 3 + yp -}— (1 + vjaT(x,9; 2) ¢ 
y 


BNp] = a & ‘ ~) Nr 
ia call B \Ox Oy Do 


b 
& (UY 5 = é,(x,¥v) (60) 


rey 1 b b 
> = 5) MNO 5 + ¢ (%,9; re = é,(x,y) (61) 
Yy 2 g 2 


where é,(x,y) is the average of the x-strain at the top and the bottom surfaces of the analog wing at the point 


(x,y) and é,(x,y) is the average of the y-strain. 


E(x,y; 2) fc 


Thus, together with Eq. (58), the above expressions give 


o(X,y; 2) = : * (ey + vi,) — 7p (Aes + vAz) — (1+ aT (xy; 5) (62) 
pb 


1—» le 


It may be proper to point out again that the quantity z in the above equation and in the subsequent equations 


is the distance from the median surface, positive if the point is below the median surface. 


one has 
E (x,y; 2) 


oy (x,y; 2) = = 
—a 


If the difference Ay,, of the shear strain and the aver- 
age shear strain y,, are defined as 


b A x 
AYry(X,¥) _ Vey X,Y; et — Tos x,y; 5 (64) 


ty. b . b 
Vr(X,y) = 9 Vay 259; ms + Fry | x,y; rs 


the shear stress in the original heated wing is calculated 


(65) 


as 
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F Fic. 3. Loads and boundary conditions of the analog test 
wing. (Numbers in parentheses correspond to numbered equa- 
tions in the text where the appropriate relation with quantities of 
the original heated wing can be found.) 





Similar to Eq. (62), 


E (é, + vez) — == (Az, + vAz,) — (1 + v) aT (x,y; 2)| (63) 
B a 


_ E(x,y; 2) & 
~ 21+ r) LB 


Tey; V3 B) 


Yru(X,¥) “ait = sia(e9) | 
ub 
(66) 


Eqs. (62), (63), and (66) allow the calculation of the 
stresses in the heated wing from test data of the un- 
heated wing. They also show the advantage having the 
constants (8/a) and uw smaller than unity, so that the 
strains of the unheated wing are magnified in computing 
the stresses in the heated wing. Then, for a specified 
load on the heated wing, the test load on the unheated 
wing is reduced and the analog wing will not be over- 
strained. Fig. 3 gives a summary of the loading of the 
unheated test wing, together with the equation numbers 
of relations that specify the quantities in terms of the 
quantities of the original heated wing. 


ALTERNATE TEST PROCEDURE FOR THIN SOLID WINGS 


In the previous section, the simulated tests on the 
isothermal wing are proposed with simultaneous load- 
ing both in the plane of the median surface and in lateral 
This results in a rather complicated load 
To simplify 


pressure. 
system at the free boundary of the wing. 
the load system at the boundary, the loading in the 
plane of the median surface and the lateral loading can 
be separated—that is, two separate tests, properly 
correlated, are made on the analog wing, and the 
stresses in the heated wing are obtained by a synthesis 
of the test results. 
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To determine first the sectional forces N,, N,, and 
N,,, the unheated wing is dimensioned to satisfy Eq. 
(42). It is then loaded with body forces specified by 
Eq. (57), but no lateral load is applied. At the “free 
edge’ of the wing, it is loaded according to Eq. (50). 
Then the magnitude of the sectional forces for the 
heated wing can be computed from Eq. (45) if N,, N,, 
and Ney are determined. These sectional forces can be 
calculated from the measured strains on the wing. Al- 
though theoretically there should be no bending of the 
wing, actually bending will be present because of devia- 
tions of the real wing from the idealized wing assumed 
inthe theory. The effects of such spurious bending can 
be eliminated by taking the average of the strains of the 
top and the bottom surfaces. Let the superscript 1 de- 
note the quantities induced by this loading in the plane 
of the median surface. Introduce the following nota- 


tions for the average strains: 


STRESSING 


HEATED WINGS 9 





LT T ) 
b b 
z,!1 (x,y) = - g,()) (x -3) + z,(0) (x 5) | | 
ys 2 2 | 
if b b 
(xy) = 5]? (x9 —S) + a (e955 
- 1 es I) b 
Yeu? (%Y) = 5] Ve? (HIF —5) AF Yeu? (95 5 
(67) 


Then the sectional forces V,, N,, and N,, for the heated 


wing can be calculated as 


N,(x,y) = = Dole? + vz,?] — No 
B 
. a . , 
N,(x,y) = B Dole + ve] — Nr (68) 
; al-—vp — 
N ry — 7 Do You 


When the average strains of Eq. (67) are measured, the loading in the plane of the median surface, the body 


forces, and the boundary forces can be removed. 
Let Eqs. (46) and (47) be satisfied. 
loading p on the test wing must be specified as 


pb = — wV?Mr + Aup + AN, 


The second step of the test is to load the unheated wing laterally. 
Then, in order that @ and w satisfy the same differential equation, the lateral 


ow , ww , ow 
caries ZN ey — (69) 
Ox" Oxoy Oy" 


On the “free edges’ of the wing, there is a bending moment equal to A\uM,r and an upward support force 


—hu(O0M,7/On) per unit length of the boundary. 
I : 


Let the superscript 2 denote the quantities induced by bending; then, 


. 5 : 5 rs 2a } 
Az (x,y) = a? (x,y; — _— 2) ~ OI) Ses 


b 
Ay:," (x,y) = Vxy°"" (x, o 9 i 


= (2 b ;, ) <7 { 
dy‘ | £9; > = 0(x,y) dy? (70) 


Vey” (x,y; ;) = 2b(x,y) = 
2 OxOoy 


With these strains differences between the top and the bottom surfaces, the loading Eq. (69) can be written as 


P(x,y) = — wV2?My + Aup(x,y) +d |W. ——— i ; + 


= (2) Aj; (2) : Ae (2) ™ 
b . Ny | (71) 


where the sectional forces N,, N,, and N,, are computed by using Eq. (68) and the strain data obtained during the 


first part of the test. 
summarizes the load systems for the two successive tests. 


When the strain averages of Eq. (67) and the strain 
differences of Eq. (70) are measured by the two succes- 
sive tests on the isothermal analog wing, the stresses in 
the original heated wing can be computed by using 
equations similar to Eqs. (62), (63), and (65), replacing 
y¥ and Az, Ay by Aé™, Ay™ of Eqs. (67) 
However, there is one im- 


é,¥ by &, 
and (70), respectively. 


This equation again demonstrates the feedback character of the lateral loading. 


Fig. 4 





portant difference between the test procedure of the 
present section and the test procedure of the previous 
a can be chosen independent of 8, because now 8 
is not necessarily equal to }. But a@ and X are con- 
nected because of the condition of Eq. (56). Out of 
the four constants a, 8, A, and yp, three can be chosen 
This greater degree of freedom and simpler 


one: 


arbitrarily. 
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Fic. 4. Loads and boundary conditions of the analog test 


wing, alternate procedure. (Numbers in parentheses correspond 
to numbered equations in the text where the appropriate relation 
with quantities of the original heated wing can be found. ) 


load system in each stage of tests than the test pro- 
cedure described in the previous section may prove 
to be advantageous in actual application. 


Box WING 


In the previous sections, the similarity laws for thin 
solid wings are formulated for obtaining the stresses in 
a heated wing by testing an analog unheated wing. A 
majority of actual wing structures is, however, much 
more complicated than a solid wing. For example, 
the main load-carrying member of the wing may be a 
box structure with top and bottom panels, ribs, and 
deep beams joining the top and bottom panels. Since 
only the top and bottom panels of such a box wing are 
exposed to the air stream, the aerodynamic heating of 
the structure is limited to these panels. 

To analyze such a box wing, the different structure 
elements can be broken down, and each individual part, 
such as the top panel, can be considered as a plate. 
To each individual part then, the theory developed in 
the earlier sections of this paper can be applied. The 
thickness 6 for the analog unheated plate for each part 
of the structure can be specified by the method de- 
scribed previously. The only restriction is that the 
four constants a, B, \, and uw must be the same for all 
structure elements. The loads on the individual un- 
heated plates corresponding to the different structure 
elements can be determined according to equations 
given in the earlier sections. At the boundary or the 
junction of the different unheated analog panels, there 
One type of load comes 


are now two types of load. 
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: 
from the similarity law. 
previously as the load system for ‘‘free edges.’’ The 
other type of load comes from the junction condition 
that the deformations of the different parts of the 
structure must fit together. If the analog panels are 
put together, then the junction loads will be automati- 
Therefore the pro- 


This load system is specified 


cally supplied by the structure. 
cedure to use the method of similarity is as follows: 

The first step is to determine the thickness 5 of the 
analog isothermal box wing by Eqs. (54) and (56). Ina 
box structure, there is appreciable influence on the 
sectional forces in the median surface of the plates 
because of bending deflection. Therefore the separation 
of bending from the “‘extensional loads”’ in the plane of 
the plates is no longer appropriate. The complete load 
system for the analog wing must be applied simultane- 
ously. This consists of the body force loads in the plane 
of the plates, the lateral pressure p, and the ‘“‘free 
edge’ loads at the junction of the different parts of the 
structure. The stresses in the original heated box wing 
can then be calculated by using Eqs. (12), (63), and (65) 
with the measured strains. 


DISCUSSION 


In the previous sections, it is shown that the problem 
of stressing a heated wing can be solved by tests on an 
analog unheated wing through the similarity laws. 
The temperature distribution in the wing is assumed 
to be given by a theoretical calculation using aerody- 
namic heating data and the theory of heat transfer in 
The loading system required for the 
Furthermore, the 


solid materials.! 
analog wing is indeed complicated. 
lateral pressure load p specified by either Eq. (59) or 
Eq. (71) involves the concept of feedback in that / is 
partially dependent upon the measured strains and 
cannot be predetermined. This is, perhaps, a novel 
concept in structure testing. The possible justification 
for this great complication in structural testing is the 
serious difficulty of simulating aerodynamic heating in 
a structures laboratory and of measuring strains at high 
temperatures. 

The basic principle of the theory of similarity law 
can be traced to the well-known fact that the general 
three-dimensional thermal stress problem can be re- 
duced to a problem in isothermal material by intro- 
ducing appropriate fictitious three-dimensional body 
forces and surface forces on the boundary of the body. 
The general three-dimensional formulation of the 
“analogy” is, however, not useful for structure testing 
purposes, because there is no possible method of apply- 
ing three-dimensional body forces. If one dimension 
of the body is small, such as a thin plate or a thin shell, 
then the body force is two-dimensional. The loading for 
the analog isothermal structure can be done, although 
not easily, as shown in the previous discussion. From 
argument, it is clear then that the 
for thermal stresses can 


this general 
present “similarity theory” 
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certainly be extended to any thin elastic shells, although 
the actual execution of the analog testing may be more 
difficult than the flat plate case studied here. 

To avoid the task of testing a heated wing, a test wing 
at room temperature is proposed as an analog. In a 
sense then, the whole concept of the method of simi- 
larity is that of analog-machine computing. There is, 
however, an advantage of the present method over the 
recognized machine computer in that the main physical 
member of the problem, the elastic plate, remains and 
is not replaced by an approximate system such as elec- 
tric network. The unheated analog wing is thus the 
closest analog the aerodynamically heated wing can 
have and is the most accurate analog. The price of 
this accuracy and detailed reproduction of the original 
problem is the complicated test setup required. How- 
ever, an enterprising structures test engineer will prob- 
ably welcome such challenge to his ingenuity and de- 
mand on his technical skill. 


APPENDIX—YOUNG’S MODULUS PROFILE 
By writing f(7) as 
f(n) = V (1/12) (a3/d) g(n) (72) 


Eq. (56) can be written as 


1 "74 “1 i 
in| =| f fia) dn || fn) dn| - 
| [to n , “ nn” an — 
*1 2 
| f fla)-ndn | (73) 
—4 


-~ 


Now any continuous function f(n) for —1 <7 < 1can 
be expanded into a series of Legendre polynomials 
P,,(n)*—.e., 


* The author is deeply indebted to Prof. A. Erdélyi, of the 
California Institute of Technology, for suggesting this method of 


solving the problem. 
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pm a,P,(n) (74) 


n= 0 


f(n) = 


where a,,’s are constant coefficients. Then‘ 


1 ra 
J t(n) dn / f(n)Po(n) dn =2ao 
e —1 —j} 
1 me 
’ 2 
f(n)ndn = f(n)Pi(n) dn = =a 
J foonan J_ SOP) ay = 5% 


° a 
/ f(n)n? dn = sf fin)| Pat) + ; Pata) in = 
1 : 1 - 


9) ») 

> = a2 + ao 

0 0 
Therefore Eq. (73) can be written as 


2 far l fa,\* = 
Qa? =1+ 5 ( ‘) — ( ) (75) 
» \do oOo \do 


The restriction of Eq. (73) is thus simply a relation 
among the first three coefficients do, a, and d» of the ex- 
pansion Eq. (74). The later coefficients, a,, m > 3, are 
Therefore the f(y), although not en- 
Hence, 


a 


entirely free. 
tirely arbitrary, has a wide degree of freedom. 
it is probable that the Young’s modulus profile will 
conform to the condition of Eq. (56) for the similarity 


theory up to a high order. 
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The Effect of the Torsional Rigidity of a Single 
Stiffener on the Buckling Characteristics of a 
Panel Subjected to Axial Compression with 
Large Deflections’ 
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SUMMARY 


The solution of von Karman’s fundamental equations for large 
deflections of plates is presented for the case of a rectangular 
panel (2:1) with a single stiffener along its center parallel to the 
applied edge compressive load for the conditions where the stiff 
ener has very low and then very high torsional rigidity relative to 
the plate. In both cases, the bending rigidity of the stiffener is 
taken of high order of magnitude compared with the plate itself 
and is assumed to remain unbuckled under all edge compressive 
loads considered. 

After solving for the deflection coefficients defining the deflec- 
tion patterns, it is possible to compare the critical buckling loads, 
effective plate widths, unit plate shortenings, membrane stresses, 
and bending stresses for the two cases. 

On the basis of these comparisons, predictions may be made 
concerning the increase of permissible load-carrying capacity of 
panels containing rigid stiffeners against those containing non- 


rigid stiffeners. 


NOMENCLATURE 


a = length of square side of plate 

a, = effective width of plate 

AR RR 

d,, 5, pq? 
a = coefficients in stress function 

c = moment arm for equivalent normal pressure, 
p 

D = flexural rigidity of plate = Eh®/12(1 — yp?) 

E = Young’s modulus 

F = stress function 

G = shear modulus 

h = plate thickness 

J = polar moment of inertia 

ae . = edge moment coefficients 

m, = edge bending moment per unit length about 
x-axis 

m,n, p,q,7,S = positive integers used as subscripts 

, a = pressure coefficients for normal pressure 

u,v, Ww = displacements of a point in the x, y, and z 
directions, respectively 

W = lateral displacement-thickness ratio = w/h 

ae coefficients in infinite series for w 
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%, 2 = coordinates of a point measured from origin 
at corner of plate 

tes'8y = unit plate shortening in x and y directions 

o,', ¢,', Tz,,/ = median-fiber stresses (membrane) 

o,",0,',Tz,y = extreme-fiber bending and shearing stresses 

€.', €y') Yx.y = median-fiber strains 

Go Saw extreme-fiber bending and shearing strains 

T,, Fy, = average edge compressive stresses 

Cer = critical compressive stress 

€cr = critical compressive strain 

mm = Poisson’s ratio (taken as 0.316) 

Vi _ o# 9 o% ‘ o* 


tr 2 t 
ox! Ox*0y? oy! 


INTRODUCTION 


N THE DESIGN OF THIN PLATES that bend under lateral 
I and/or edge loading, formulas based on the Kirch- 
hoff theory, which neglects stretching and shearing in 
the middle surface, are satisfactory provided that the de- 
flections are small compared with the thickness. If 
deflections are of the same order as the thickness, the 
Kirchhoff theory may yield results that are consider- 
ably in error, and a more rigorous theory that takes 
account of deformations in the middle surface should 
therefore be applied. The fundamental equations for 
the more exact theory have been derived by von Kar- 
man.! The analysis of the postbuckling characteristics 
of plates in this paper are based upon these nonlinear 
equations of von Karman. 

The numerical solutions presented herein are based on 
the trigonometric series method of attack introduced by 
Levy.’ 

The critical buckling loads, which are independent of 
the nonlinear terms, have subsequently been deter- 
mined by using only the linear deflection terms in the 


von Karman equations. 


EQUATIONS FOR THE DEFORMATION OF THIN PLATES 


The fundamental equations governing the deforma- 
tion of thin plates developed by von Karman! are given 
by Timoshenko? in the following form: 
O'w\? dO'w Ow 
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Vey = ph (= Ow OF ew > OF =) 
" pp D \Oy? Ox? Ox? Oy? ~ Oxdy OxXdy 
(2) 
where the median-fiber stresses are 
, OF , OF , O’F 
%. = oe = oo Se eS (3) 
oy? Ox” OxOy 
and the median-fiber strains are 
7 l ( -*) ) 
€r = a 
E \ay? ~ “ax? / | 
, l (= st) t (4 
€ = — 4) 
’ E \dx? dy? | 
, _ —2(1 +4) oF 
me E Oxoy 
The extreme-fiber bending and shearing stresses are 
‘ Eh | 2 4 4 
Or oS > Py 
2(1 — pw?) Lox? * “ ay? 
Eh Ee + a4 (5) 
g eelrat 9 2 Koo 
. 2(1 — pu?) Loy? ox? 


try’ = [Eh/2(1 + u)](0*w/dxdy) 


GENERAL SOLUTION OF THE PLATE EQUATIONS 


Fig. 1 represents a 2:1 rectangular plate with a single 
centered stiffener subjected to an end compressive load 
parallel to the stiffener. From the symmetry of the 
figure, each plate half may be considered as a square of 
length a. 

If the stiffener has a very low torsional rigidity, it will 
tend to twist freely under load so that the plate would 
be analogous to one with four edges simply supported. 
The high bending rigidity of the stiffener will resist any 
lateral deflection so that a node will exist at the center of 


the plate face. 
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The boundary conditions for this case are 


wy =O,a waa 0 | 


| 4 
- = 0; “ = 0 
Ox* x=0,a oy? y 0,a | 


This problem has been treated extensively by 
Levy.” 4 The pertinent results of his work will be pre- 
sented for comparison purposes in this paper. The 
method of attack will not be included here. 

If the stiffener has a very high torsional rigidity, it 
will resist twisting under load so that the plate will act 
as one with three edges simply supported and the fourth 


wW)x=0,0 = 9; 
(6) 


edge clamped rigidly. 
The boundary conditions for this case are 





ow ) 
w|, = O,a aii 0; Wy = 0,4 = 0: od = Q 
Ow O*w 
m = —D = 0 7 
, Fe "* 9 | 0, a 7) 


O*w Ow 
m — D + ) = 0 
Condition of Straight Edges: 


If the elongation of the plate z, in the x-direction and 
z, in the y-direction and if the edges of the plate are 
straight, then ~, and ~ must be independent of the 
position along the length or width at which they are 


measured, or 


f E (oF oF l (~) | F 
eos +r Hoy _ x 
= 0 LE \Ooy? F Ox? 2 \oOx 

_— f E es oF) 1 ()'] ; 
ty Jo LE\ax? ay?) 2\ay/ J 


Load on the Edges: 


(8) 


If the average stress in the plane of the plate is é, in 
the x-direction and 6é, in the y-direction, then 


= [(5),..7 Ge), -01 
a3, = — 
Oy] y=6 Ov/y=0 


(9 
(5) (5) | 
agy = ~ 
. Ox x=a Ox x = 0. 
General Solution: 
The lateral deflection is taken in Eq. (1) as 
= a _ marx . nTry 
w= >) D> wmasin sin —— (10) 
m=1m-=1 a a 


Eq. (10) can be made to satisfy all the boundary con 
ditious of Eqs. (7) if the value of the moment, m,, along 
edge y = a were known. 

To facilitate this solution, the required edge moment, 
m,, will be replaced by an auxiliary pressure distribution 
p(x, y) near the edge y = a, as shown in Fig. 2. This 
method was first suggested by Levy.° 
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If this pressure distribution is expressed by a Fourier 
series and the value of ¢ approaches zero, the auxiliary 
pressure is 


co 


24m, _ nwry 


(xy = > ——~ eae — (11) 


n 1 a~ a 


If m, is expressed by a Fourier series with k,, as an 
undetermined coefficient, then 


a” < ; _ Mmrx 
m, = > km sin (12) 
2 m 1 a 


Inserting Eq. (12) into Eq. (11) gives 


marx . ny 


P(x, vy) = > a Pm Sin sin (13) 
m=1n=1 a a 
where 
Pmin = (—1)"~'!nkm (14) 
From boundary conditions (7), Ow/Oy],-, = 0, it 


follows that 


p & jel ])" = Os w= 1,3... (15) 


n 1 


In order to satisfy Eq. (1) and the boundary condi- 
tions of Eqs. (9), F is taken as 


. ter es = ap MX nTy 
F = + + >» > Jua'COS cos 
2 2 m On 0 a a 
(16) 


Substituting Eqs. (16) and (10) simultaneously into 
Eq. (1) and equating coefficients of like trigonometric 
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terms on the two sides of the equation results in a rela- 
tion between the coefficients /,,,, and Wy». 
This relation can be simplified to 
E 
fan = }(m? + n?)? Le br .s.0.¢ WrsWp.¢ (17) 
where (1) the summation includes all products w,,,w,, 
for which m is equal to either the sum or difference of r 
and p, and for which v is equal to either the sum or dif- 
ference of s and q; and (2) the coefficient 6,,,.,.. is given 
by 
b, 


8D 2rspq + (r°q? + s*p?) (17a) 


where the sign before the parenthesis is positive if m is 
the sum of ry and p and x is the difference of s and g, or 
if m is the difference of r and p and 7 is the sum of sandg. 
It is negative otherwise. Whenr = pands = gq, b,.5.», 
is one-half the value given by Eq. (17a). 

It remains to show the necessary relation between the 
coefficients Wm, », Eq. (10); pm,» Eq. (13); and fn,» 
Eq. (17). 


TABLE 1 
Equations for the Stress Coefficients in Eq. (17) 





foor = (E/4)(—w, 1W1, 2 — Wi, 21, 3 — Awe, 1W3, 2 — Aws, oW3, 2s, 3) 
fo.2 = (E/32)(w1, 2? + 9ws, 1? — 2wy,, 121, 3 — 18ws. 1wWs. 3) 

fo. 3 = (E/36)(w, 1%, 2 + 9ws, 13, 2) 

fo.s = (E/128)(2uy, :w, 3 + wi, 2? + 18s, ; Ow. »2) 

fo.s = (E/100)(w, 2%, 3 + 9ws, ows, 

fo. = (E/288)(w, 32 + 9ws, 3?) 

fe,0 = (E/32)(wi, 1? + 4%, 2? + 9wi, 32 — 2wy,, 1wWs3, 1 — Swi, ows, 2 — 18wy. 3ws, 3 
fo,1 = (E£/100)(9wy, 1wi, 2 — Wy, 1W3, 2 + Zw, ow, 3 — Ww, ows, 1 — Iw, 23,3 — 49w,, Ws, 2) 
f = (E/64)(4w, 1wi, 3 — 12%, 1ws, 1) 

fe, 3 = (E/676)( —wy, 11, 2 + 251, 1w3, 2 + 19w), o%W3, 1) 

fo,4 = (E/400)( —wy, 11, 3 + 9w1, 13, 3 + 16w,; ows, » + 2Z5w, 3ws, 1) 
fo,s = (E/3,364)( —wi, ow, 3 + 81w, ows, 3 + 121, 3ws, ») 

fo. 6 = (E/1,600)(36w,, 3ws, 

fa,o = (E/64)(w, 13,1 + 41, ows, 2 + Ow, 
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fs,5 = (E/6,724)( —9w, ows, 3 — 49w), 3ws. 

fi.6 = (E/2,704)( —9wy, sw, 3) 

fe. 0 = (E/288)(ws, 1? + 4ws, 2? + Ows, 3?) 

fe.: = (E/5,476)(81ws, 1w3, 2 + 225ws, ow 

fe. 2 = (E/1,600)(36w5, ws, 3) 

fe. 3 = (E/900)( —ws, ws, 2) 

fo.s = (E/2,704)( —9ws, 1ws, 3) 


fe.6 = (E/3,721)( —9ws, ows, 3) 
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2) 


Substituting Eqs. (10), (13), and (17) simultaneously into the lateral equilibrium Eq. (2) and equating co- 
efficients of like trigonometric terms on the two sides of the equation results in the relation given by 


” 


a’ a? 


m* n* m1 
Pm.n = Din, ot ( + ) — 8hW,, . — 


a~ 


n> hr* ; 
pe T/NWm, n e + \ 4 ,» "een TrisWp,¢@ ( 1S) 
a~ a 


where (1) the summation includes all products /,, .w,,, for which m is equal to either the sum or difference of r 
and p, and for which v is equal to either the sum or difference of s and g; and (2) the coefficients d,, ,, », , are given 


by 


ie S,p.q@ 


and are twice this value if either 7 or s 1s zero. 


The first sign is positive if either r — p = morifs —q 


in all other cases. 


+(rq + sp)? ifr ¥ 0,5 #0 (18a) 


The signs in Eq. (18a) are determined by the following rules: 
= n but not if both conditions are true. 


It is negative 


The second sign is positive if m is the sum of r and p and 1 is the difference of s and q; or if m is the difference 


of rand p and n is the sum of s and q. 


It is negative otherwise. 


The boundary condition for straight edges, Eqs. (8), can be shown to be satisfied by the chosen expcessions for 
F and w by substituting Eqs. (10) and (16) into Eqs. (8) and performing the indicated differentiation and integra- 


tion. 
amount 


orl G,a 


‘e=-- sae i 


I E 


while the edges y = 0 and y = a separate by an amount 


G,a 4 G,a 
ee . Mo} 
E E 


After simplification, requiring the use of Eq. (17), it is found that the edges x = 0 and x = a separate by an 


m2 is 
> Dd mn, 2" (19a) 
Sd m=12=1 
ae . 
> Dd 82m. 2° (19b) 


od m=12=1 


Since the right-hand sides of Eqs. (19a) and (19b) are independent of x and y, the boundary conditions for straight 


edges are satisfied. 


Equations (15), (17), and (18) each represent infinite sets of equations. 
A good approximate solution can, however, be obtained by assuming that only 


should satisfy all of these equations. 


An exact solution for this problem 


a limited number of w,,, , coefficients differ from zero and that, of these, only a few will have large enough values 


so that they cannot be neglected where third-order terms appear. 
The values of f», , in Eq. (17) for these coefficients are listed in Table 1. 


Eq. (18) represents a doubly infinite family of equa- 
tions. In each of the equations of the family, the co- 
efficients f,,, may be replaced by their values listed in 
Table 1. 

The resulting equations will involve the unknown nor- 
mal pressure coefficients, p,, », the cubes of the deflec- 
tion coefficients, w,,,, and the known average edge 
compressive stresses, ¢, and é,. For this problem 4, is 
equal to zero. The number of these equations is equal 
to the number of unknown deflection coefficients w,,. , 
and is consequently six. 

By utilizing Eq. (14), the six unknown pressure co- 
efficients, Pm, », can be replaced by two moment co- 
efficients, k; and ks. 

These six equations are solved for the linear-deflec- 
tion term in terms of the remaining terms and the values 
used for the first three terms in the summations of the 
two of Eqs. (15) obtained by taking m = 1 and m = 3 
(the remaining deflection coefficients being taken as 
zero). The two cubic equations listed in Table 3 result. 
Thus, there are eight simultaneous cubic equations with 


eight unknowns. A direct solution is virtually impossi- 
ble so that iterative methods must be employed. 

It will be noted that the equations in Tables 2 and 3 
are cubics, and their solution therefore gives three values 
for each of the deflection coefficients, w,,.,. Some of 
these values correspond to stable equilibrium, whereas 
the remaining values are either imaginary or correspond 
to unstable equilibrium. Fortunately, if the equations 
in Tables 2 and 3 are solved by a method of successive 
approximations, the successive approximations ap- 
proach values corresponding to stable equilibrium. 

The final solution for the deflection coefficients and 
edge moment coefficients are shown in Table 4. 

With the values of the deflection coefficients known, 
the horizontal unit plate shortening may be calculated 
from Eq. (19a). The variation of the horizontal unit 
plate shortening with the average edge compressive 
stress is represented in Fig. 3. 

The effective width of plate, defined by 


Ae = [(82/ Ocr.z)/(Er/€cr.z) Ja (20) 
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Fic. 3. Variation of horizontal unit plate shortening with aver- 
age edge compressive stress. 
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Fic. 4. Variation of effective width ratio with horizontal unit 
plate shortening. 


is plotted against horizontal unit plate shortening in 
Fig. 4. 

The error resulting from the omission of more deflec- 
tion terms in the solution can be judged to some extent 
by comparison of the critical stress constant ¢,a°/Eh? = 
6.45, Table 4, with the value 5.24 of the same constant 
computed from the tables in reference 6 for a square 
plate having one edge parallel to the load clamped and 
the other edges simply supported. Since the results 
given in reference 6 correspond to an exact solution of 
the differential equations in the small deflection range 
that governs buckling, the difference is a measure of the 
error in the present solution. It is estimated that the 
error in the large deflection range is slightly larger than 


TABLE 4 

G,a*/Eh? Wi, 1 Wi, 2 Wi, 3 W3, 1 Ws, 2 

6.450 0 0 0 0 0 
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9.451 1.5 0.4919 —0.1721 0.0866 0.0351 
11.526 2.0 0.6312 —(). 2457 0.1588 0.0693 
14.162 2.5 0.7554 —0.3293 0.2407 0.1175 
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The outstanding advantage resulting from the use of 
a stiffener of high torsional rigidity with respect to the 
plate as illustrated in Fig. 3 is the substantial increase 
in the critical stress constant, ¢,a°/Eh?, from the value 
of 3.66 to 6.45. 
per cent in plate area is possible for constant thickness 


In this manner, a reduction of about 33 


plates having a stiffener of high torsional rigidity in 
preference to plates having a stiffener of low torsional 
rigidity. 


Figs. 3 and 4 show a reduction of unit plate shortening 
of about 25 per cent and an increase of effective plate 
width of about 20 per cent for the plate containing a 
torsionally rigid stiffener. 


The quantitative value of the results is limited by the 
number of terms used in the determination of the de- 
flection patterns. With the use of six deflection co- 
efficients, a reasonably close approximation of the de- 
flection patterns is attained. In both cases, the w;, ; co- 
efficients are far greater numerically than any of the 
other coefficients. 
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The nonlinear character of plate problems makes 
possible, under certain conditions, more than one stable 
solution. In this paper, only one real root for each de- 
flection coefficient resulted from the solution of the 
cubic equations. However, if a lateral pressure existed 
in the conjunction with the edge load, a load might be 
reached where one stable deflection configuration ‘‘oil- 
cans’’ into another stable configuration. This stability 


problem might well be treated ia some future work. 
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Boundary-Laver Stability and Transition 
in Subsonic and Supersonic Flow 
A Review of Available Information with 


New Data in the Supersonic Range* 


CARL GAZLEY, JR.7 


General Electric Company 


SUMMARY 


Accurate prediction of the aerodynamic characteristics of a 
body at high flight speeds requires knowledge of the characteris 
tics of the existing boundary layer. Although the effects of com 
pressibility on the characteristics of both laminar and turbulent 
boundary layers are known rather well, little information has been 
available as to how the transition between the two types of flow 
is affected by compressibility and other factors. Even in subsonic 
flow, the effects of such variables as free-stream turbulence, sur 
face curvature, pressure gradient, surface roughness, and surface 
temperature are known only qualitatively. 

While the theory of laminar boundary-layer stability yields 
the conditions necessary for instability, it does not permit pre 
diction of the transition point. This theory is, however, useful in 
indicating possible effects of the several variables on transition 
It has recently been extended to compressible flows with heat 
transfer and with pressure gradient 

This paper presents the available data in subsonic and super 
sonic flow for the effects of free-stream turbulence, surface curva 
ture, pressure gradient, surface roughness, surface temperature, 
and Mach Number on the transition position. New supersonic 
data are included from rocket flights, firing-range tests, and wind 
tunnel tests. The effects of the several variables on transition are 
compared with the theoretically predicted effects on boundary 
layer stability 

In general, the experimental data for transition confirm the 
trends indicated by the stability theory. Perhaps the most sig 
nificant deviation from the trend expected from the stability 
theory is the relative insensitivity of supersonic transition to 
surface-temperature variation. The transition Reynolds Number 
range appears to increase with increasing Mach Number in both 
firing-range tests, where relatively low surface temperatures 
occur, and in wind-tunnel tests, where relatively high surface 


temperatures occur 


NOMENCLATURE 


B = width of roughness element in flow direction 
k = height of roughness element above surface 
f, = function of 


function of 


/ = characteristic dimension of turbulence-producing de 
vice 

MI = Mach Number 

p = pressure 

R = Reynolds Number 
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7 = absolute temperature 
u = velocity in x-direction (time average) 
u’ = fluctuating component of u 
ux = friction velocity, ./7,/p 
x = distance in flow direction from leading edge or nose 
apex 
6 = boundary-layer momentum thickness 
\ = Pohlhausen (or Weil) parameter 
uw = absolute viscosity of fluid 
vy = kinematic viscosity of fluid 
p = fluid density 
T = shearing stress 
Subscripts 
aw = adiabatic-wall condition 
? = roughness element position 
t = transition position 
zw = wall condition 
Cs = based on x 
i] = based on 6 
| = condition just outside boundary layer 
©» = free-stream condition 


INTRODUCTION 


— POSITION OF BOUNDARY-LAYER TRANSITION is 
one of the factors necessary for the prediction of 
aerodynamic performance in high-speed flight. This is 
apparent when one considers that the condition of the 
boundary layer (i.e., whether laminar, turbulent, or 
separated) determines the skin-friction drag, the base- 
pressure drag, the form drag of the afterbody, and the 
effectiveness of control surfaces that are partially im- 
mersed in the body boundary layer. The effects of the 
boundary layer on the various types of “‘pressure’’ drags 
and control-surface effectiveness become more im- 
portant as the flight speed increases. This is due to the 
fact that boundary-layer thickness increases with Mach 
Number (at constant Reynolds Number). 

When the present investigation was initiated about 
2 years ago, there was little quantitative information 
available regarding the effects of various quantities on 
boundary-layer transition, especially in supersonic flow. 
Researches in subsonic flow had shown that the transi- 
tion Reynolds Number is affected by surface tem- 
perature, surface roughness, pressure gradient, surface 
curvature, and free-stream turbulence. Investigations 
in transonic flow indicated that shock-wave-boundary- 


layer interaction also influences transition. However, 
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only the qualitative effects of these variables were 
known, little quantitative information was available 
about the effect of these variables, either singly or in- 
terrelated. Quantitative information about boundary- 
layer stability resulted from the theories of Tollmien! 
and Schlichting? and was confirmed by the experiments 
of Schubauer and Skramstad.* However, this theory 
and its extensions only yield the Reynolds Number 
above which disturbances can be amplified in the 
boundary layer and does not give the Reynolds Number 
where the disturbances have been sufficiently amplified 
to cause turbulence. The theories are practically use- 
ful, however, for indicating possible qualitative effects 
of various quantities on transition. 

Because of the many factors that affect boundary- 
layer transition, the experimental study of the effect of 
one of these variables on transition is difficult. For ex- 
ample, the effects of free-stream turbulence in wind- 
tunnel studies tend to mask out the effects of other 
variables. This is especially true of supersonic flow, 
since little is known about turbulence in supersonic flow 
or its effects; furthermore, the turbulence characteris- 
tics of most supersonic tunnels have not been deter- 
mined. The desirable method of experiment is, of 
course, tests in flight; however, they are expensive and 
can be used only in cases where other methods are not 
possible. An attractive possibility appears to be the use 
of firing-range tests for transition study. There is no 
free-stream turbulence present, and it is possible to 
evaluate the transition position from spark shadow- 
graphs of the projectile in flight. 

This paper presents the available data in subsonic 
and supersonic flow for the effects of free-stream tur- 
bulence, surface curvature, pressure gradient, surface 
roughness, surface temperature, and Mach Number. 
New supersonic data are included from rocket flights, 
firing-range tests, and wind-tunnel tests. The effects of 
the several variables on transition are compared with 
the theoretically predicted effects on boundary-layer 
stability. 


THEORY OF LAMINAR BOUNDARY-LAYER STABILITY 


The theory of the stability of laminar flows was first 
formulated mathematically by Rayleigh.‘ Successful 
solutions for the laminar boundary layer on a flat plate 
were obtained by Tollmien' and by Schlichting.? The 
results of this theory indicate that instability in the 
laminar boundary layer originates in small disturbances 
that are either damped or amplified depending on their 
frequency and the Reynolds Number. This selective 
amplification is predicted only if the Reynolds Number 
is greater than a certain value, called the ‘minimum 
critical Reynolds Number.’ This Reynolds Number is, 
somewhat smaller than the transition 


” 


of course, 


Reynolds Number, since the regular disturbances in the 
laminar boundary layer must be amplified to a magni- 
tude sufficient to produce the irregular oscillations 
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characteristic of turbulence. For several years, no ex. 
perimental evidence was found of the laminar boundary- 
layer oscillations predicted by the Tollmien-Schlichting 
theory. 

In 1935, a somewhat different theoretical approach 
was suggested by Taylor.° The theory of Taylor, 
originally thought to be in opposition to that of Toll- 
mien and Schlichting, involves the supposition that 
turbulence in the free stream outside the boundary 
layer imposes local adverse pressure gradients on the 
laminar boundary layer. These local adverse gradients 
are thought to cause local boundary-layer separation 
with subsequent transition. 

The apparent difference between these two theories 
was finally resolved by the experimental studies of 
Schubauer and Skramstad,* at the National Bureau of 
Standards. In connection with a study of wind-tunnel 
turbulence, these investigators observed the regular 
laminar boundary-layer oscillations predicted by Toll- 
mien and Schlichting in a wind tunnel having an ex- 
ceptionally low turbulence level. The predicted ranges 
of damping and amplification were closely confirmed. 
It was found that the oscillations and their damping or 
amplification to eventual turbulence could be detected 
only when the magnitude of the free-stream turbulent 
fluctuations is less than about 0.001 of the mean 
velocity. For larger magnitudes of the free-stream tur- 
bulence, the boundary-layer oscillations are difficult to 
identify because of the near coincidence of their appear- 
ance and the point of transition to turbulence. 

With the aid of the Schubauer-Skramstad experi- 
ments, the Tollmien-Schlichting theory and the Taylor 
theory can be reconciled, and a unified picture of the 
mechanism leading to transition in the boundary layer 
can be obtained. In the case of low free-stream 
turbulence (u’/u, < 0.001), any small disturbances 
that are initially present may be due to either internal 
or external disturbances such as surface roughness, 
vibration, noise, or free-stream turbulence. The ampli- 
fied oscillations cannot be classified as turbulence be- 
cause of their regularity; it is only when these oscilla- 
tions or waves become large enough to roll up into eddies 
that turbulence is initiated. The eddies, formed close 
to the solid surface, break away and dissipate in an 
irregular turbulent motion. The process is well illus- 
trated in the interferometer photographs of Eckert® and 
is somewhat analogous to the growth and breaking of 
waves on an air-water interface. In case of a high 
free-stream turbulence level (u’/u,; > 0.001), the free- 
stream fluctuations are large enough to produce local 
adverse pressure gradients that cause local separation 
with subsequent transition in the manner described by 
Taylor's theory. 

Since the Tollmien-Schlichting oscillations have been 
confirmed experimentally, revisions and extensions of 
the theory have appeared. Lin’ has made some re- 
visions and mathematical clarifications. Schlichting 
and Ulrich’ and Hahneman, Freeman, and Finston® ex- 
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BOUNDARY-LAYER 


tended the theory to flows with pressure gradient and 
found that a positive pressure gradient destabilized, and 
a negative pressure gradient stabilized, the boundary 
layer. The minimum critical Reynolds Number is 
shown in Fig. 3 as a function of the pressure-gradient 
parameter. Lees and Lin"™, !! and Van Driest!? have 
studied flat-plate compressible flows with heat trans- 
fer. Weil'® studied the effect of pressure gradient in 
compressible flow. 

While the theories of the Tollmien-Schlichting type 
permit prediction of the effect of various factors on 
laminar boundary-layer stability and prediction of the 
Reynolds Number below which disturbances are not 
amplified, they give no indication of the transition Reyn- 
olds Number. Presumably, the original small disturb- 
ances must be amplified to the point where they be- 
come unstable. While the theory permits calculation of 
the initial amplification rate, it is applicable only when 
the disturbances are infinitely small. That this ampli- 
fication takes considerable time and distance is ex- 
emplified by a comparison of the minimum critical 
Reynolds Number with the transition Reynolds Num- 
ber. For incompressible flow along a flat plate, the 
minimum critical Reynolds number is calculated as 
60,000; this figute was confirmed by Schubauer and 
Skramstad,* who also found the Reynolds Number 
range of the transition region to be 2,800,000 to 3,900,- 
000. For this case, then, there is about a sixty-fold 
difference between the Reynolds Numbers at which 
amplification begins and at which turbulence begins. 


Liepmann'‘ has made an interesting attempt to com- 
pute the transition Reynolds Number using the maxi- 
mum initial rates of amplification given by the Toll- 
mien-Schlichting theory. He estimates that transition 
to turbulence will begin when the apparent shear stress 
of the amplified disturbances becomes equal to existing 
laminar shear stress. The resulting relationship gives 
the transition Reynolds Number as a function of the 
magnitude of the initial disturbance having the fre- 
quency of maximum amplification. However, this 
latter quantity is dependent on the magnitude and fre- 
quency of the free-stream turbulence, of the surface 
roughness, and/or of noise and vibration; the quanti- 
tative dependence of this factor on these variables is 
not known. Therefore, Liepmann’s calculation cannot 
be compared directly with experiment; however, for 
reasonable values of the initial disturbances, the pre- 
dicted ratio of transition Reynolds Number to minimum 
critical Reynolds Number is of the same order as that 


found experimentally. 


Emmons" has recently evolved a theory of transition 
based on the idea that every point on a body is turbu- 
lent part of the time; the fraction of the time that a 
given spot is turbulent is derived from probability con- 
siderations. The results of this theory, while not ena- 
bling prediction of the absolute value of transition range 
of Reynolds Numbers, give the variation of skin friction 
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or heat-transfer rate in the transition region and the 
effect of flow geometry (flat plate, cone, swept wing) on 
the transition. 


EXPERIMENTAL STUDIES OF THE EFFECTS OF VARIOUS 
FACTORS ON TRANSITION 


Because of the importance of the position of transi- 
tion to the practical aerodynamicist, many experi- 
mental determinations have been made of the effects of 
various factors on transition, and many experimental 
techniques have been devised to detect transition. 
Most of the available experimental data are in the 
low subsonic-flow range but are useful for giving 
an insight into the effects of such factors as pressure 
gradient, free-stream turbulence, surface roughness, 
etc. While these data enable qualitative knowledge of 
the effect of a given variable upon transition, they do not 
usually allow a quantitative correlation. This is gen- 
erally due to the fact that other variables, also affecting 
transition, were not controlled or measured during the 
experiments. Thus it is frequently necessary to com- 
pare data taken under different conditions, although in 
recent years the appreciation of the sensitivity of transi- 
tion to a great variety of factors has yielded several 
valuable studies in which extraneous factors have been 
carefully controlled. Notable among these are the in- 
vestigations of Liepmann, at the California Institute of 
Technology,'* ' and Dryden and coworkers, at the 
National Bureau of Standards.* These investigators 
carefully controlled the free-stream turbulence, pres- 
sure gradient, and surface condition and used a hot- 
wire anemometer to detect transition. This is probably 
the most sensitive technique used for detecting the prog- 
ress of transition. Other techniques involve the dif- 
ference in velocity distribution, skin friction, transfer 
rates, or temperature recovery factor in the laminar and 
turbulent boundary layers. 


Free-Stream Turbulence-—As noted above, Taylor*® 
has developed a theory for the effect of free-stream 
turbulence on transition. The result of this theory in- 
dicated that the transition Reynolds Number is a func- 
tion of the intensity and scale of the free-stream tur- 


bulence: 
1 
Ru = fi ((u’/m)(x/l) | (1) 

where 

u’ = free-stream velocity fluctuation 

u, = free-stream velocity 

x = distance downstream from turbulence-pro- 

ducing device 
1 = characteristic dimension of turbulence-produc- 


ing device 


This relationship was confirmed experimentally for 
transition on spheres at the National Bureau of 
Standards.” Recent studies of turbulence characteris- 
tics'® have indicated that the intensity and scale of iso- 
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Fic. 1. The effect of free-stream turbulence on boundary-layer 
transition in subsonic flow. 


tropic turbulence can be related approximately as 

x/l =~ u/u’ (2) 
Thus Eq. (1) may be simplified to 

Ru = fo (u'/m) (3) 


Of course, these relationships would be expected to hold 
only under conditions where the scale of the turbulence 
is of the same order of magnitude as the boundary- 
layer thickness. For example, atmospheric turbulence 
is presumably of such a large scale that it would have 
only a negligible effect on transition. 

The available data for the change in transition 
Reynolds Number with free-stream turbulence in- 
tensity are shown in Fig. 1 for the case of subsonic flow 
over smooth flat plates aligned with the stream. It is 
seen that the general form of Eq. (3) is confirmed and 
that the variation of free-stream turbulence can 
change the transition Reynolds Number by as much as 
thirtyfold. It is interesting to note that even the higher 
levels of free-stream turbulence (fluctuations of about 
3 per cent of the mean velocity) do not produce transi- 
tion Reynolds Numbers below the minimum critical 
ralue calculated by Lin.’ 

Surface Curvature.—The effect of surface curvature on 
flow stability has been analyzed by many investigators. 
Qualitative considerations by Rayleigh* and Prandtl!® 
indicate that flow along a concave wall is less stable 
than flow along a convex wall. Quantitative studies by 
Schlichting” and Géortler?! of two-dimensional dis- 
turbances indicate the opposite. Recent studies by 


Taylor’? and Gértler** of three-dimensional disturb- 
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ances indicate that flow along a concave wall is less 
stable than flow along flat or convex walls. 

The effects of surface curvature on boundary-layer 
transition observed by Liepmann' confirm approxi- 
mately the theory of Goértler**—1.e., that the boundary- 
layer stability is dependent on the value of the param- 
eter RoV 6/r which, for zero pressure gradient, is pro- 
portional to R,“V x/r. The experimental transition 
Reynolds Numbers observed by Liepmann, for zero 
pressure gradient, are shown in Fig. 2 in comparison 
with lines of constant Goértler parameter. The data for 
low turbulence levels (u’/u, = 0.0006, 0.002) give a 
Gortler parameter of about 9.0 and for a higher tur- 
bulence level a value of about 5.5. The effects of curva- 
ture in the region of convex curvature seem to be negli- 
gible, as Goértler’s theory indicates. 

Pressure Gradient.—The effects of pressure gradient 
on laminar boundary-layer stability have been com- 
puted by Schlichting and Ulrich® and Hahneman, et al.,’ 
for incompressible flow and by Weil!* for compressible 
flow. The results of these calculations, shown in Fig. 
3,7 indicate that a negative pressure gradient tends to 
stabilize the boundary layer while a positive gradient 
has the opposite effect. The predicted effect of pressure 
gradient is greatest at a Mach Number of 1.5 and is al- 


Tt These results are shown in terms of the momentum-thickness 
Reynolds Number, since the x Reynolds Number is without 
significance in flows with pressure gradient. 
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Fic. 2. The effect of surface curvature on boundary-layer 
transition in subsonic flow. Data of Liepmann" for dp/dx = 0 
and u’/u, = 0.0006. 
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most negligible at a Mach Number of 4. General ex- 
perience in subsonic flow confirms qualitatively the 
effects predicted theoretically—e.g., transition fre- 
quently occurs shortly after the minimum pressure 
point on an airfoil. Furthermore, the experiments of 
Schubauer and Skramstad? confirmed qualitatively the 
predicted effects of pressure gradient on the damping or 
amplification of laminar boundary-layer oscillations. 
However, only one set of systematic experiments for 
the effects of pressure gradient on transition is known 
(those of Liepmann"™), and none is known in com- 
pressible flow. The results of Liepmann, also shown in 
Fig. 3, follow the same trend as that predicted by the 
theory. As is to be expected, the transition measured 
by Liepmann occurs at higher Reynolds Numbers than 
the theoretical minimum critical Reynolds Numbers. 

Surface Roughness.—Both the surface finish and any 
single surface discontinuities are known to affect transi- 
tion in the boundary layer. Small irregularities in the 
surface finish presumably contribute to the initiation of 
oscillations in the laminar boundary layer which are 
selectively amplified until transition occurs. How- 
ever, no known quantitative information exists on the 
relationship between the scale and characteristics of a 
surface finish and the resulting boundary-layer oscilla- 
tions. Larger single roughness elements may bring 
about an almost immediate transition if they are large 
enough to produce a turbulent wake. Considerable 
data exist for the effects of two-dimensional roughness 
elements on transition, such as wires and ridges, and 
a few data exist for three-dimensional elements, such as 
rivet heads. Most of the data lies in the subsonic- 
flow region, although a few data are becoming available 
in supersonic flow. 

The available data for the effects of two-dimensional 
roughness elements cover a variety of shapes and sizes 
mounted on bodies of revolution, flat plates, and air- 
foils. Measurements have been made in a water 
tunnel,?4 several wind tunnels,” and in flight.” Fage and 
Preston*! formulated an empirical criterion for the size 
of wire necessary to cause immediate transition—i.e., 
transition directly behind the wire. Tani, et al.,* pro- 
posed a somewhat similar correlation for the effect of a 
wire. Goldstein” proposed a similar criterion for the size 
of a surface wire which has no effect on the transition 
position. Fage* tested various types of two-dimensional 
elements and devised another empirical correlation for 
the effect of the roughness elements on the transition 
position. All of these correlations can be reduced to the 


form 
Ru = f(k?/Bx) (4) 
where 
R,, = transition Reynolds Number 
k = height of roughness element 
B = width of roughness element 


x, = distance to transition point from leading edge 
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Fic. 3. The effect of pressure gradient on boundary-layer 
stability and transition. Transition data of Liepmann! for sub- 
sonic flow; stability calculations of Hahneman, et al.,’ for in- 
compressible flow and of Weil'® for compressible flow over 5 per 
cent thick biconvex airfoil. 


The available data for the effect of two-dimensional 
roughness elements on transition are shown in Fig. 4, 
together with the several proposed correlations. It is 
apparent that Fage and Preston’s correlation for wires** 
fits the data well in the lower Reynolds Number range, 
while Fage’s correlation for a variety of other shapes” 
follows the data at higher Reynolds Numbers. It is 
apparent that any general correlation should asymp- 
totically approach the smooth-surface value of the 
transition Reynolds Number for extremely small rough- 
ness elements. None of the heretofore proposed correla- 
tions has this characteristic. A proposed correlation, 
covering the entire Reynolds Number range and ap- 
proaching the smooth-surface transition Reynolds 
Number at small roughness values, is shown by the solid 
lines in Fig. 4. Admittedly, the degree of correlation is 
not too good but is perhaps as good as can be expected 
considering the variety of conditions under which the 
data were procured—i.e., with different body geome- 
tries, pressure gradients, and degrees of free-stream tur- 
bulence. Undoubtedly the effect of roughness is de- 
pendent on the free-stream turbulence. Most of the 
high Reynolds Number data shown in Fig. 4 was taken 
either in flight or in low turbulence-level tunnels; how- 
ever, it may be seen from Fig. 1 that smooth-surface 
transition may occur at lower Reynolds Numbers if the 
free-stream turbulence is appreciable. Therefore, a 
rariety of curves might be drawn for the roughness 
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Fic. 4. The effect of single spanwise roughness a, on 


boundary-layer transition in subsonic flow. 

Line A-A, limit of buffer layer in resulting turbulent boundary 
layer, Ceriiadtaee »y = 30. 

Line B-B, suggested correlation of Fage and Preston,** ku, /v 
= 20. 

Line C-C, suggested correlation of Tani, et al.,%* ku, /vy = 13. 

Line D- D, suggested correlation of Goldstein, 7 ku, /v = 7. 

Line E-E, suggested correlation of Fage.* 


Symbol shapes denote roughness shapes: (©) = wire;?4~7* — 
= smooth bulge or half wire; 77 ~ = smooth hollow;> — = 
flat ridge.™ Solid points indicate x, = x; flagged points in- 


dicate flight data. 


correlation for various turbulence levels. Two are 
shown in Fig. 4, one for low turbulence and one for 
“average”’ turbulence. 

The surprising fact deducible from this correlation is 
that the relative position of the roughness element does 
not affect the position of boundary-layer transition. 
Fage” indicates that such an effect does exist only when 
the roughness element is located at relatively small 
values of x. However, the inclusion of this factor into 
the correlation, as suggested by Fage, does not appre- 
ciably improve it. The apparent fact that transition is 
not appreciably affected by the relative location of the 
roughness element indicates that, as the roughness ele- 
ment is moved downstream, the decreasing disturbance 
due to the element and the increasing inherent insta- 
bility of the laminar boundary layer are combined to 
produce a constant instability due to the roughness 


element. It may be noted that the critical size of rough- 


ness elements indicated in Fig. 4 is quite a different 
function of the Reynolds Number than is the bound- 
ary-layer thickness, the critical roughness size varying 
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from the same order of magnitude to a small fraction of 
the boundary-layer thickness. 

Data on the effects of three-dimensional roughness 
elements are rather meager; however, the data of 
Loftin® and of Williams and Brown’ for rivets agree 
reasonably well with the correlation for two-dimen- 
sional elements when the rivet-head diameter is used 
for the factor B. A small amount of data for back- 
ward lapped joints” indicates that the joint has ap- 
proximately the same effect as a wire whose diameter is 
one-half the joint height. The data of Holstein® for 
various grades of emery cloth on a flat plate correlate 
reasonably well in the form shown in Fig. 4 but lie 
somewhat below the line shown for two-dimensional 
elements, approaching Goldstein’s line at low Reyn- 
olds Numbers. Unfortunately, Holstein’s 
appear to be affected by a variation of tunnel turbu- 


results 


lence level occurring during a change in tunnel pressure 
level. 

There is little information available about the effects 
of roughness in supersonic flow; a few data are shown 
in Fig. 5 and tend to confirm the trend indicated in 
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Fic. 5. The available data for the effect of roughness elements 
on boundary-layer transition in supersonic flow. Symbol shapes 
denote source of data: A = Station H of V-2 No. 27; © = 
cone in Aberdeen Firing Range; [_] missile model in wind tunnel. 
Shading of points indicate position of transition: open points 
denote no transition, shaded points denote transition down- 
stream of roughness (x; > x,), and solid points denote transition 
at roughness (x» = x,). Arrows on open and solid points in- 
dicate probable position of a corresponding shaded point. 
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Fig. 4 with the Mach Number effect shown in Fig. 8 
superimposed. 

Surface Temperature.—The theory of boundary-layer 
stability'®~'* indicates that surface cooling has a stabi- 
lizing effect while surface heating has a destabilizing 
effect. The results of Van Driest’s recent calcula- 
tions!® are shown in Fig. 6. It is apparent that the 
predicted effect of surface temperature is appreciable, 
although not so great as originally calculated by Lees." 
The available experimental data*!~** for the effect of 
surface temperature on transition are also shown in 
Fig. 6. While the observed effect of surface tempera- 
ture on the transition Reynolds Number is seen to be 
qualitatively similar to the predicted effect on the 
minimum critical Reynolds Number, the quantitative 
effect is somewhat less. The data shown are internally 
consistent as far as the effects of surface temperature 
are concerned; however, the various sets of data are 
probably not comparable, since the turbulence level in 
the supersonic tunnels and its effect on transition are 
not known. 

Compresstbility—The theoretical works of Lin and 
Lees’ '! and of Van Driest!? for flat-plate flow indicate 
that the effects of both Mach Number and heat trans- 
fer on laminar boundary-layer stability are appreciable, 
as is indicated in Fig. 6. Van Driest’s results (which 
are extensions and corrections to the earlier work of 
Lees) indicate that, between Mach Numbers of | and 9, 
the boundary layer may be completely stabilized by 
means of a considerable degree of cooling. Fig. 7 shows 
Van Driest’s computed minimum critical Reynolds 
Number as a function of Mach Number and surface 
temperature. If the surface temperature is maintained 
in constant ratio to the free-stream temperatures, an 
increase in Mach Number is predicted to increase the 
stability of the boundary layer; if the surface tempera- 
ture is maintained in constant ratio to the adiabatic- 
surface temperature, an increase in Mach Number is 
predicted to destabilize generally the boundary layer, 
except for the intermediate Mach Number range when 
VT < 3Tqw/4, where complete stability is predicted to 
exist. As noted above (Fig. 6), the meager transition 
data available confirm the predicted trend of the effect 
of surface temperature on stability, although the transi- 
tion trend is not so pronounced as the predicted sta- 
bility trend. However, no systematic experimental 
investigations have been made of the combined effects 
of surface temperature and Mach Number. 

Because of the effects of free-stream turbulence on 
transition (Fig. 1), the study of the effects of Mach 
Number on transition in wind tunnels is of questionable 
value unless the free-stream turbulence is suppressed 
to a point where it no longer influences transition 
Unfortunately, the turbulence characteristics of super- 
sonic wind tunnels are not known; in fact the charac- 
terization of turbulence in compressible flow has no 
firm foundation. Furthermore, the effects of free- 
stream turbulence on transition have not been estab- 
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Fic. 6. The effect of surface temperature on boundary-layer 
stability and transition on nominally smooth surfaces. The 
transition data shown at the top of the figure indicate the be- 
ginning or middle of the transition region, except for the data of 
Higgins and Pappas,*? which show the beginning and end of the 
transition region. The stability calculations'* are for a flat 
plate. 
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Fic. 7. The effect of Mach Number and surface temperature 
on boundary-layer stability. The calculated stability trends! 
are for constant values of 7./Taw (solid lines) and for constant 
values of 7/7; (dotted lines). The numbers adjacent to the 
curves denote values of these temperature ratios. 
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Fic. 8. The effect of Mach Number on boundary-layer tran- 

sition in firing range (7. ™ 7;) and in rocket flight (7. ™ 1.3 7). 

= cones in Aberdeen Firing Range; [] = V-2 No. 27;% 

A V-2 No. 21; VY V-2 No. 12; and © rocket flight. Open points 

denote beginning of transition region; solid points denote the 
end. 
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Fic. 9. The effect of Mach Number on boundary-layer transi- 
tion in wind tunnels. <1 = flat plates and airfoils;5* *% %, 1 33, 
3, 41, 42 () = bodies of revolution*! 4% #4 (including unpublished 
G-E data). Open points denote beginning of transition, shaded 
points denote average transition point, and solid points denote 
end of transition. 


lished in supersonic flow. Because of these considera- 
tions, measurements in flight or in a firing range ap- 
peared to be attractive possibilities. Some flight skin- 
temperature data had already been obtained by the 
author’s organization on V-2 rockets;® these data, 
when used to compute convective heat-transfer rates, 
indicated the transition point of the boundary layer. 
Some additional skin-temperature data later became 
available from other rocket flights. In the rocket 
flights, the measured skin temperature lay between the 
temperature outside the boundary layer and the adia- 
batic-wall temperature. Most of the measurements 
were made on the smooth conical noses of these 
rockets; thus the effects of pressure gradient, as well as 
free-stream turbulence, were eliminated from these re- 


sults. Additional free-flight data were obtained from 
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firing-range shadowgraphs of cone-cylinder models 
loaned by A. C. Charters and J. D. Nicolaides, of the 
Ballistics Research Laboratory, Aberdeen Proving 
Ground. By careful examination of these shadow- 
graphs it was possible to detect the apparent position of 
the beginning and end of transition on the conical nose 
of the models, where again a zero pressure gradient 
existed. In these firings, the time of flight was so short 
that the surface temperature remained close to ambient 
temperature. Effects of surface finish, however, may 
have been present in the firing-range tests, since the 
models were only machine-finished. The transition 
Reynolds Number ranges estimated from the rocket 
flights and the firing-range tests are shown in Fig. 8 asa 
function of the local Mach Number outside the bound- 
ary layer. Although the data show considerable 
scatter, there appears to be an increase in stability with 
increasing Mach Number. This is the trend predicted 
for these surface temperature conditions by the sta- 
bility theory, although the effect of Mach Number on 
transition is not so great as the predicted effect on 
stability. 


Although wind-tunnel transition data are suspect 
because of the unknown effects of tunnel turbulence 
level, the available data are shown in Fig. 9. These 
data were obtained from continuous-flow tunnels where 
the models had presumably reached equilibrium tem- 
perature. For this surface-temperature condition, the 
stability theory predicts decreasing stability with 
Mach Number. The wind-tunnel data, however, 
show a trend similar to that shown by the flight and 
firing-range data. Thus it may be tentatively con- 
cluded that the trend of boundary-layer transition with 
surface temperature is considerably less than the pre- 
dicted trend of stability. This agrees with the rather 
slight trend with surface temperature shown in Fig. 6. 
Both sets of data (Figs. 8 and 9) indicate an increase in 
transition Reynolds Number with Mach Number. 
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Fic. 10. Comparison of laminar boundary-layer stability 
theory with boundary-layer condition on smooth nose cone of V-2 
rocket in flight. Data of Fischer and Norris® for Stations A, C, 


and G on the smooth nose cone of V-2 No. 27. 
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BOUNDARY-LAYER 


These must be considered as tentative conclusions, 
however, because of the possible effects of surface 
roughness on the firing-range results and of tunnel 
turbulence in the wind-tunnel results. The only data 
that may be considered free from these restrictions are 
those of Fischer and Norris*® for Stations A, C, and G 
on the nose cone of V-2 No. 27. These data are shown 
in Fig. 10 compared with the predicted region of com- 
plete stability calculated by Van Driest.'? It is seen 
that stabilization of the boundary layer takes place 
close to the theoretical prediction. This is not true, 
however, for stations preceded by surface roughness, as 
Fig. 11 indicates. Station H was located on the nose 
cone opposite Station G but was preceded by an !/s-in. 
trip wire. Stations K and M were located back on the 
ogival nose of the missile in a region of negative pres- 
sure gradient but were preceded by door joints and other 
The wind-tunnel and firing- 
it is seen that 


surface discontinuities. 
range data are also indicated in Fig. 11; 
the appearance of turbulence in the wind-tunnel tests is 
compatible with the theory, since they lie in the “un- 
stable’ region. However, in the firing-range tests, 
turbulence appears in the “‘stable”’ region; this may be 
due to the machine-finish of the firing-range projectiles. 


CONCLUSIONS 


A review of the available boundary-layer transition 
data and a comparison with theory have yielded certain 
conclusions as to the effects of several variables on 
transition. These conclusions are summarized quali- 
tatively in Fig. 12. The predicted effects of several of 
these variables on laminar boundary-layer stability are 
shown, as well as their experimental effect on transition. 
It is apparent that, in general, the experimental transi- 
tion trends follow the predicted stability trends. The 
only exception appears to be the effects of surface 
temperature and Mach Number on supersonic transi- 
tion, which are much less than, and sometimes opposite 
to, the predicted effects on stability. 

It is evident that, while the qualitative effects of the 
several variables on transition are well established in 
subsonic flow and partially in supersonic flow, quantita- 
tive information is rarely available. This is particu- 
larly true of the combined effects of several of the vari- 
ables that affect transition. Additional data on the 
effects of surface roughness, surface temperature, and 
pressure gradient in supersonic flow would be particu- 
larly valuable. The use of a firing range appears 
especially attractive for this work, at least until the 
turbulence level in supersonic wind tunnels and its 
effects on transition are established. 
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Equilibrium Temperature of an Unheated 
Icing Surface as a Function of Air Speed’ 


BERNARD L. MESSINGERt 
Lockheed Aircraft Corporation 


ABSTRACT 


The thermal analysis of a heated surface in icing conditions has 
been extensively treated in the literature... Except for the work 
of Tribus,! however, little has been done on the analysis of an 
unheated icing surface. This latter analysis is significant in the 
design of cyclic thermal deicing systenis that are attractive for 
small high-speed aircraft for which thermal anti-icing require- 
ments have become severe. 

In this paper, a complete analysis of the temperature of an 
unheated surface in icing conditions is presented for the several 
significant régimes (i.e., less than 32°F., at 32°F., and above 
32°F.) as a function of air speed, altitude, ambient temperature, 
and liquid water content 

The results are presented in graphical form and permit the 
rapid determination of surface temperature for a wide range of 
Curves are presented to determine the speeds beyond 
Curves are also presented to 


variables. 
which no ice accretion will occur. 
indicate the surface temperature and the rate of ice sublimation 
which takes place when an ice-covered surface emerges into clear 
air. 

One significant result of this study is the introduction of a new 
basic variable referred to as the “freezing-fraction,’’ which de- 
notes the proportion of the impinging liquid which freezes in the 
impingement region. The fact that some of the liquid does not 
freeze in the impingement region tends to explain the observed 
variation in ice formation shape with temperature, speed, and 
water catch. 

New test data obtained at Mt. Washington, N.H., for stagna- 
tion-point surface temperatures of an unheated plastic cylinder in 
natural and artificial icing conditions are included in the Ap- 
pendix. These data substantiate the validity of the assumptions 
made in the theoretical analysis. 


(1) INTRODUCTION 


ne MAY BE PROTECTED against icing in two 
basic ways. The ice may be entirely prevented 
from forming (anti-icing) as, for example, by thermal or 
chemical means; or it may be allowed to form to a 
tolerable thickness and periodically removed (deicing) 
as in the case of the inflatable boot or the cyclic applica- 
tion of heat. Until recently, almost all aircraft em- 
ployed either an inflatable boot type of deicing or a 
thermal system of anti-icing, the latter being the pre- 
ferred method. 


Presented at the Design and Structures Session, Annual Sum- 
mer Meeting, I.A.S., Los Angeles, June 27-28, 1951. Revised 
and received July 2, 1952. 

* The author is indebted to M. M. Friedlander for assistance 
in preparation of the curves and to P. J. Valentine for suggestions 
relative to the “freezing fraction”’ concept. Both are members of 
the Thermodynamics Group at Lockheed Aircraft Corporation. 


t Department Manager, Thermodynamics. 
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The advent of transonic aircraft speeds and very thin 
airfoils has made the application of thermal ice preven- 
tion extremely difficult or impractical, and, as a result, 
increased attention has been given to the possibilities 
offered by the cyclic thermal deicing method. The 
publication of reference 1 has been an important factor 
in directing attention to the significant economy of 
energy that is possible by the use of cyclic thermal de- 
icing as compared to the continuous thermal anti-icing 
method. A more recent report, refererce 2, contains 
experimental results of preliminary N.A.C.A. research 
with cyclic electrothermal deicing using a relatively low- 
intensity system. By using a short duration, high- 
intensity ‘‘on’’ period, the average energy can be of the 
order of one-tenth to one-twentieth of energy required 
for continuous anti-icing. This large reduction in en- 
ergy requirements makes electrical power appear to be 
an ideal source of heat for a cyclic thermal system. 
These general characteristics have been experimentally 
verified by numerous flight and laboratory tests re- 
cently conducted by the Lockheed Aircraft Corpora- 
tion. 

Of the thermal types, the continuous anti-icing sys- 
tem has predominated in the last decade. Hence, the 
theoretical studies and methods of analysis developed 
in this period have been almost entirely concerned with 
the energy balance that exists when an anti-iced surface 
is maintained at a temperature considerably above the 
freezing point, generally at about 100°F. Little if any 
attention (except in reference 1) has been directed to- 
ward determining the temperature of an unheated icing 
surface. 

If, in designing a cyclic thermal deicing system, high- 
intensity heating (about 30 to 40 watts per sq.in.) is 
used, a short ‘‘on’’ period as low as 2 sec. may be em- 
ployed. This permits a ratio of “‘off’’ to ‘‘on”’ time of 
as much as 50 or 75 to 1. With this large a ratio of 
“‘off’’ time, the inttial temperature of the iced surface 
just prior to an ‘‘on”’ period is essentially the equilibrium 
temperature of an unheated airfoil surface. Although 
the exact mechanics of the ice removal process in a 
cyclic thermal system is not fully understood, it is be- 
lieved that this initial surface temperature has a signifi- 
cant bearing on the net energy required. Other factors 
that may be involved are related to the effect of ice 
thickness and liquid surface tension on the removal 


process. It was pointed out in reference 1 that one 
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advantage of the deicing method is that a high rate of 
water catch results in a higher initial temperature in the 
ice impingement area due to the release of the latent 
heat of fusion, whereas an anti-icing system is severely 
taxed at high rates of water impingement because not 
only is the heat of fusion unavailable but, in addition, 
the entire latent heat of evaporation must be supplied, 
as well as the usual loss by convection heat transfer. 
By using high-intensity thermal deicing, it is possible to 
achieve nearly instantaneous, dry removal of ice, thus 
avoiding the effect of ‘‘runback’”’ and consequent re- 
freezing. 

During the early phases of the development of air- 
craft ice-protection systems, it was usually assumed 
that, when the day came that airplanes would fly fast 
enough to result in stagnation temperatures of 32°F. 
or more, the problem of icing would be eliminated. 
As that day approached, however, it became increas- 
ingly obvious that the surface temperature rise in icing 
conditions would be different from that which would 
prevailin dry air. It is now evident that ice can be de- 
posited on a surface even though the speed corresponds 
to a clear-air surface temperature of more than 32°F. 

One of the early attempts to analyze the effect of 
free moisture on the kinetic heating of the boundary 
layer was made by J. K. Hardy, a visiting British anti- 
icing specialist who was temporarily stationed at the 
Ames Aeronautical Laboratory of the N.A.C.A. The 
results of his work were reported in 1944 in reference 3 
and somewhat later in references 4 and 5. In reference 
4, Hardy introduced the concept of a ‘‘datum tempera- 
ture,’ which he considered to be a basic reference param- 
eter for use in computing the unit heat- and mass- 
transfer rates for a heated surface in icing flight. This 
“datum temperature” parameter was adopted and used 
in many subsequent reports and articles on ice protec- 
tion.°-* In some of these articles (particularly refer- 
ence 6), the suggestion is made that the “‘datum temper- 
ature’ cannot be measured in an icing cloud because of 
the release of the heat of fusion. In the opinion of the 
author of the present paper, there are certain fallacies 
in the manner in which the “‘datum temperature” param- 
eter has been employed, and these will be discussed in 
detail later in the article. 

It is the purpose of this paper to present a complete 
analysis of the conditions that govern the equilibrium 
temperature of an insulated, unheated surface exposed 
toicing. The analysis is, in part, based on the approach 
begun in reference 1 but goes beyond the scope of that 
work. The latter considered icing conditions for a sur- 
face temperature above and below 32°F. but not at 
32°F. The present analysis includes an important 
range of conditions at a 32°F. surface temperature using 
the concept of the ‘freezing fraction.” In extending 
the results of this analysis to the problem of airfoil ice 
protection, some consideration must be given to the 
limitations of some of the simplifying assumptions in- 
volved in certain of the parameters that are employed. 
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In that range of conditions in which a surface temper- 
ature of more than 32°F. is attained during icing, an 
analysis has also been developed using a parameter 
that will be referred to as the ‘evaporation fraction,” 
This latter analysis includes conditions up to and in- 
cluding the speed for complete evaporation of all inter- 
cepted water. 

Also of interest to the aircraft designer are those fac- 
tors that govern the ability of the airplane to rid itself 
of ice accretions collected on unprotected surfaces 
after it has emerged from icing conditions into clear air. 
For this reason, an analysis is presented which permits 
the determination of (a) the rate of sublimation of ice 
from a surface at less than 32°F. in clear air, (b) the 
speed required to attain a surface temperature of 32°F, 
on an ice-covered surface in clear air, and (c) the rate 
of melting from a wet surface of ice at 32°F. in clear 
air. 


(2) SyMBoLs 


A = surface area, sq.ft. 

B = ambient absolute pressure, in. of mercury 

b = “relative heat factor,’’ Rwcw/fe, dimensionless 

Cs = unit heat capacity of ice, B.t.u./lb. °F. 

by = unit heat capacity of air, B.t.u./lb. °F. 

ti = unit heat capacity of water, B.t.u./lb. °F. 

6m ‘= linear melting rate, in. per hour 

5s = linear sublimation rate, in. per hour 

En = water droplet catch efficiency, dimensionless 

Te = unit convection conductance, B.t.u./hour ft.? °F. 

g = gravitational constant, 32.2 ft./sec.? or lbs. per slug 

vi = density of ice, lbs. per cu.ft. 

J = mechanical equivalent of heat, 778 ft.lbs. per B.t.u. 

L = latent heat of vaporization of H2O, B.t.u. per Ib. 

Ly; = latent heat of fusion of HO, B.t.u. per Ib. 

L, = latent heat of sublimation of ice, B.t.u. per Ib. 

m = evaporation fraction, dimensionless 

n = freezing fraction, dimensionless 

Fr = vapor pressure of atmospheric moisture, in. of mercury 

P,; = vapor pressure over ice at ¢s,, in. of mercury 

P.. += vapor pressure over water at é,,, in. of mercury 

q = rate of heat flow, B.t.u. per hour [subscripts defined 
in paragraphs (3.1), (3.2), and (6.2)] 

r = “recovery” factor applying to kinetic heating, dimen- 
sionless 

Ry» = unit rate of water catch, lbs./hour ft.? 

Rw, = unit rate of evaporation of water, lbs./hour ft.? 

Rom = unit rate of melting of water, lbs./hour ft.? 

t= ambient free-stream temperature, °F. 

ts. = equilibrium surface temperature, °F. 

0 = heat balance parameters, °F. [defined in paragraph 
(4)] 

V. = free-stream velocity, ft. per sec. 


(3) Moprs oF ENERGY TRANSFER TO AN UNHEATED 
IcING SURFACE 


(3.1) Equilibrium Surface Temperature Less Than 32°F. 


For an icing surface having a steady-state tempera- 
ture, f:,, of less than 32°F., the actual value of ¢,, will 
be largely dependent on the result of an energy balance 
involving the simultaneous interchange of the following 
quantities (see Fig. 1): 
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ANALYSIS OF AN UNHEATED 


(3.1.1) Heat lost by convection: 
Ye = fA (ts — ta) 
(3.1.2) Heat Jost by sublimation: 
gs = 2.90L,f.A[(Ps, — P.)/B]* 
(3.1.3) Heat lost due to being absorbed by the 


warming of the impinging subcooled liquid: 
Qu = RyAc,(32 — t,) 


(3.1.3) Heat gained due to the release of the latent 
heat of fusion as the subcooled liquid impinges and 
changes to the solid state at 32°F. and, then, as a solid, 


cools to ¢;,: 
7 = RyA[144 + c,(32 — ¢,,)] 


(3.1.5) Heat gained due to the viscous or frictional 
heating in the boundary layer: 


Wo = fA(rV,,?/2gJcp) Tt 


(3.1.6) Heat gained equivalent to the kinetic energy 
of the liquid particles as they strike the icing surface: 


gx = RyA(V,*/2gJ) 


Because, in this case, the value of ¢,, is assumed to be 
less than 32°F., all of the impinging water is considered 
to have solidified so that no mass transfer by direct 
“blowoff”’ of liquid takes place and the accretion rate 
will be R,,. In addition, the effect of heat loss by radia- 
tion has been entirely disregarded, since it is relatively 
insignificant. 


(3.2) t,, Equal to 32°F. 


For a surface having an equilibrium temperature of 
32°F., the energy balance is basically similar to that 
outlined above, except that an additional variable (here- 
after referred to as the ‘“‘freezing fraction,’’ ”) is in- 
volved in energy quantity (3.1.4). When operating 
conditions result in the value of ¢,, just reaching 32°F., 
n will be equal to 1.0, which signifies that a// of the im- 
pinging liquid is converted to ice. When conditions 
correspond to the release of a somewhat greater net 
amount of heat to the ice formation (as with increased 
speed), part of the impinging water will not solidify but 
will persist as a 32°F. liquid, so that the value of m, the 
“freezing fraction,’ will be /ess than 1.0. Whenever n 
is less than 1.0, the heat quantity (3.1.2) corresponds 
to the latent heat of evaporation rather than heat of sub- 


* The constant 2.90 is an empirical factor relating mass transfer 
to convection heat transfer. In reference 1 this constant appears 
erroneously as being equal to 2.29 because of a numerical error 
in the derivation. Also, to be completely accurate, the values of 
P., and B should be replaced by their corresponding local values 
just outside the boundary layer. Thus the areas of low static 
pressure promote increased evaporation and high static pressure 
areas suppress evaporation. The magnitude of this effect is con- 


sidered negligible in this paper. 
t This term is usually combined with that of paragraph (3.1.1) 
to give the met convection loss. 
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Fic. 1. Modes of energy transfer for an unheated airfoil in 
icing conditions. 


limation, even though a large part of the water may be 


present as ice. 
The energy balance for a 32°F. surface temperature, 


t,., in icing, therefore involves: 


(3.2.1) Heat lost by convection: 
Ge = f,A(32 — t,) 
(3.2.2) Heat lost by evaporation (unless m = 1): 


de = 2.90L. fA [(Pso — P..)/B) 


Heat lost due to warming of the liquid: 


Qo = RyAc,(32 — t,) 


(3.2.4) Heat gained due to the latent heat of fusion: 
dr = 144nR,A 
(3.2.5) Heat gained due to viscous or frictional heat- 
ing: 


de» = fA(rV,2/2¢Icp) 


(3.2.6) Heat gained due to the kinetic energy of the 
impinging liquid: 


ge = Ry,A(V,,?/22J) 


As the value of m becomes smaller, it is probable that 
a significant part of the impinging liquid will leave the 
surface by “‘blowoff”’ or ‘‘runoff.”” When the value of 
n becomes zero, none of the impinging water freezes 
and the surface may be at 32°F. or above. Note that 
the actual ice accretion rate in this range will be equal 
to nRy. 


(3.3) t., Greater Than 32°F. 


For above-freezing surface temperatures, at least 
some of the impinging water will probably run off the 
trailing edge of the airfoil during actual icing conditions. 
At the same time, however, part of the impinging water 
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will leave the airfoil by evaporation, and that fraction 
that is evaporated in the impingement area is denoted 
by the symbol, m, and can be determined analytically 
as will be shown in a later section. At certain condi- 
tions corresponding to m = 1.0, the net energy gained 
by the surface will permit complete evaporation of the 
impinging liquid; hence, there will be no runback and 
presumably the surface will approach a dry condi- 
tion. 

The energy balance for a surface temperature, ¢,,, 
greater than 32°F. in icing, therefore involves: 

(3.3.1) Heat lost by convection: 


Ye = fA(ts — te) 
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(3.3.2) Heat lost by evaporation: 
2.90L.f-A[(P:, — P.,)/B] = 
Heat Jost due to warming of the liquid: 


Iw = RyA Cuts, = t..) 


mR,L,A 


= 


(3.3.3) 


(3.3.4) Since none of the liquid freezes, there is no 
term representing the heat of fusion. 


(3.3.5) Heat gained due to frictional heating: 
g = f.A(rV.,?/2gJcp) 
(3.3.6) Heat gained due to the kinetic energy of the 
impinging liquid: 


gx = R,A(V,,2/2¢J) 


(4) EVALUATION OF THE BASIC HEAT BALANCE EQUATIONS 


(4.1) When t,. < 32°F. 


Arranging the terms of paragraph (3.1) in a single heat balance equation yields: 
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tse for 20,000-ft. altitude. 


v 
2g] 


The last term can be evaluated at c, = 1.0 and c; = 
0.47, 


R,,A | et? — 4.) — 144 — of{38 — t,.) - 


... + R,A[0.47t,, — t. 29J)] 


— 127 — (V,.? 


By introducing* the dimensionless ratio b = RyCi/fe 
and three new groupings of the variables—namely, 


6; = t,,(1 + 0.47b) + (2.90L,P,/B) 
6. = t,(1 + b) + (2.90L,P,,/B) + 1276 
03 = [(r/cp) + b](V.2/22J) 


the combined equation then becomes 
gq = foA(@, — A — 83) 


For the case of an unheated surface under equilibrium 
conditions, g = 0 and therefore 6; = 42 + 43. 

For any given altitude, the function @; can be plotted 
in terms of ¢, and b; 62, in terms of ¢,, and b; and 6; for 
an arbitrary value of r [r = 0.875, a compromise value 
for laminar (0.85) and turbulent (0.90) flow], in terms 
of V, and 6. Such plots are shown in Figs. 2, 3, and 4. 
These 6 functions have the dimensions of temperature 
°F. and, except where they are independent of altitude, 
have been plotted for 20,000 ft. 


* These parameters were first presented in reference 1. 
+ A more precise form for r is [1 — (V2/V,2)(1 — Pr”)], 
where n = !/, for a laminar and m = !/; for a turbulent boundary 


layer. 
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(4.2) When t,, = 32°F.and0 <n <1 


For this case, the terms of paragraph (3.2) can be con- 
veniently correlated by means of a second group of 
6’ functions: 


6,’ = 32(1 + b) + (560/B) 
(based on LZ, = 1,075 for ¢,, = 32°F.), 
t.(1 + b) + (3,120P,,/B) + 144nb 
03 = [(r/cp) + 6](V.*/2gJ) 
and also, in this case, since g = 0, 


6,’ = 62’ + 6,’ 


6,’ 
63’ 


II 


For any given altitude 6,’ is a function only of } as 
shown in Fig. 5, and 6,’ can be plotted in two steps as 
indicated in Figs. 6 and 7. 


(4.3) When t,, > 32°F. and n = 0 


The terms of paragraph (3.3) can be expressed by a 
third group of @” functions: 


0,” = t,,(1 + 6) + (3,100P,/B) 
(based on L, = 1,070 when 32°F. < t,, < 50°F.), 
6.” = t,(1 + b) + (3,100P,,/B) 
(based on L, = 1,070 for 32°F. < t,, < 50°F.), and 
6,” = 0,’ = 03 = [(r/cp) + b](V2/2g¢/) 
and again 


A; ” =n 0.” ao 6,” 
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As indicated in Fig. 8, both 6,” and @,” can be repre- 


sented by a single plot. 


(4.4) When t,, > 32°F. and n < 1 


In order to determine the extent of the evaporation 
in the impingement area when ¢,, > 32°F. and n < 1, 
the following energy balance can be solved in terms of 
m, the evaporation fraction: 

The energy required for evaporation: 


Je = mMR,~L.A 
The energy available for evaporation: 


ge’ = 2.9f-L-A[(Ps,. — P..)/B] 
Equating the two forms, 
- a s) 
8 B 


2.9 ./P.—- P. 

a 
Ry B 

This latter equation can also be used to determine 

the speed beyond which a surface exposed to icing con- 

ditions will be maintained ice-free and dry (m = 1): 


P,, = b(B/2.9) + P, 


After finding P,,,, the corresponding value of ¢,, can 
be found in any table of steam properties. The speed 
corresponding to this value of ¢,, can then be deter- 
mined by means of Fig. 4 from @;” = 6,” — @,”. 
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(5) EQUILIBRIUM SURFACE TEMPERATURE IN ICING 
CONDITIONS AS A FUNCTION OF THE SIGNIFICANT 
PARAMETERS 


(5.1) Characteristics of the Parameters 

By means of the graphical presentations of the @ 
functions, it is a simple matter to evaluate the heat 
balance equation and determine the equilibrium tem- 
perature ¢,, for any part of the range of flight speeds in 
icing conditions regardless of whether this temperature 
is above, at, or below the freezing point, provided, 
however, that the local value of b can be determined 
with reasonable accuracy. 

The dimensionless ratio 6 (which is equal to Ry.c./f,) 
can be considered a measure of the ratio of the sensible 
heat-absorbing capacity of the impinging water per 
unit of surface area to the unit convective heat-dissi- 
pating capacity of the same surface. It has been pro- 
posed that this ratio be called the ‘relative heat fac- 
tor.” 

No general statement can be made about the varia- 
tions of the value of 6 except that, for a given position 
in the impingement area of an airfoil and for a constant 
air speed, it will vary directly as the liquid water con- 
tent. 

Another characteristic of the parameter, 0, is that it 
tends to increase with altitude, since R,, increases and 


f, becomes smaller for a given speed. 


The variation of 6 with air speed depends, in turn, on 
the relationship between the variation of R, and f, with 
air speed for a given position on the airfoil. For an 
area in which the initial catch efficiency is close to 
100 per cent, the value of R, will vary directly with 
speed; if the initial catch efficiency is low, R, will in- 
crease at a more rapid rate than the air speed. If the 
external boundary is locally laminar, f, will vary as the 
square root of the air speed; if turbulent, it will vary 
as the 0.8 power of the air speed; in the transition re- 
gion, it may vary as an air-speed power function greater 
than 1.0. It is therefore apparent that no simple single 
relationship will express the variation of 6 with air 
speed. 

By making a survey of the local values of R, and f, 
around the periphery of a typical thin airfoil, it is pos- 
sible to investigate the range of variation of } in the 
impingement area. Based on a laminar conductance, 
the maximum value of } is about three to four times the 
stagnation point value and seems to occur at about 5 
per cent chord on the lower surface for an airfoil angle 
of attack at 4°. A typical stagnation point value for } 
is 0.5 based on a 6-ft. chord N.A.C.A. series 65-210 
airfoil at 20,000 ft., 435 knots, 0°F., 30 micron drop 
size, and a liquid water content of 0.5 Gm. per cu.m. 

A study of the @ function characteristics for equilib- 
rium conditions yields the following relationships: 

(5.1.1) 6, cannot be greater than 6, since 6; cannot 
be negative. 
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Fic. 9. tse versus V for b. = 0.1 and b = 0.2 at 10,000-ft. 
altitude. 


(5.1.2) When 6: = 6, a purely hypothetical condi- 
tion prevails—.e., 63 = 0(V. = 0)—and yet both water 
impingement and heat transfer are assumed to exist. 
This relationship, however, does establish the minimum 
t,, that can exist at subfreezing conditions for a given 
value of b. [See paragraph (5.2) and Fig. 10a for a de- 
tailed analysis of an actual condition. | 

(5.1.3) If 6. is greater than that value of @,, which 
corresponds to 32°F., than m must be less than 1.0— 
i.e., @ wet icing surface exists—and the @’ functions 
should be used. 

(5.1.4) Ifm = Oandt,, = 32°F., for a given value 
of /,,, 92” can be read directly and then subtracted from 
6,” (corresponding to ¢,, = 32°F.) to give the required 
values of @;” and thus the speed for a surface that is wet 
but on which no ice is depositing. This method was 
used to construct Fig. 11 from Fig. 4 and Fig. 8 for 
20,000 ft. and from the equations for 0,” and 6,” of para- 
graph (4.3) for 10,000 and 30,000 ft. 

(5.1.5) If the speed falls in the range between the 
two conditions defined by (5.1.3) and (5.1.4) above 
(i.e., 0 << < 1), the following steps can be used to de- 
termine the freezing fraction, m, for given values of 
V., t., 6, and altitude or to determine V for given 
values of n, ¢,,, 6, and altitude: 

(5.1.5.1) For the given value of altitude and ), de- 
termine the corresponding value of 6,’ using Fig. 5. 

(5.1.5.2) If V., is given, determine @;’ using Fig. 4. 

(5.1.5.3) Calculate 6.’ = 6,’ — 63’. 
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For given value of ¢,, determine 6.’ — 


(5.1.5.4) 
(144nb) using Fig. 6, and then calculate (144d) using 
§,' as obtained in (5.1.5.3). 

Using Fig. 7, determine m for the given 


(5.1.5.5) 
value of 0. 

(5.1.5.6) If m is a known value and it is necessary 
to determine the corresponding value of V «, then step 
(5.1.5.2) should consist of determining 6.’ by adding 
(144nb) obtained from Fig. 7 to the value of 6.’ — 
(144nb) obtained from Fig. 6. 

(5.1.5.7) Using Fig. 4, V can then be determined 
from 6;’, which is the sum of 6,’ and 64’. 


(5.2) Surface Temperature Characteristics 


In order to illustrate the general trend of ¢,, as a 
function of airplane altitude, speed, water catch rate, 
and surface conductance, Figs. 9 and 10 have been con- 
structed by use of the @ functions of Parts (4.1), (4.2), 
and (4.3). 

Although each of these plots of ¢,, versus air speed is 
shown for constant values of }, it should be noted that b 
will vary somewhat with speed as indicated above. It 
has been suggested by J. P. Lewis, one of the coauthors 
of references 2, 6, and 9, that, since Figs. 9, 10, and 12 
do not present a realistic picture of the true variation 
of surface temperature with air speed, it would be de- 
sirable to construct these curves for constant values of 
liquid water content as would be the case in actual icing 
flight. Unfortunately, it is only possible to construct 
such a plot for a given location on a given configuration 
of aerodynamic body. In order to illustrate the differ- 
ences between the general plot using b as the independ- 
ent parameter and a specific plot using liquid water 
content as the independent parameter, Fig. 10a has 
been prepared based on the stagnation region of a 3-in. 
diameter cylinder. The lower part of the figure shows 
the variation of b with air speed for this point on the 
cylinder, and the shaded area of the upper part of the 
figure indicates the nature of the discrepancy that re- 
sults from the use of the simplifying assumption that 5 
remains constant. As noted on this figure, the accu- 
racy of the catch efficiency in the higher speed range is 
doubtful, but it can be seen that a large variation in b 
has little influence on /,, at the higher speeds. 

A line of constant liquid water content (= 0.4 Gm. 
per cu.m.) has also been added to the lower part of Fig. 
12 to illustrate the same effect at sea level. 

Fig. 10 shows that the effect of a high rate of water 
catch (b = 0.5) is to maintain the impingement area at 
32°F. over a wide range of speeds, in this case from 185 
to 587 knots. Note that, while a temperature of 32°F. 
is attained at only 185 knots, there is at least some ice 
accretion all the way up to 587 knots. Even at 395 
knots for b = 0.50, Fig. 10 indicates that 60 per cent of 
the deposit is ice. There is as yet no experimental in- 
dication that, under these conditions, the liquid water 
present (40 per cent at 395 knots) acts in any way that 
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would assist aerodynamic forces in blowing away the 
ice formation. In fact, a few high-speed icing flights 
have been unofficially reported for at least moderate 
icing conditions, at which time heavy, rapidly forming 
deposits of ice were observed. The liquid water present 
when 2 < 1 may form runback ice beyond the impinge- 
ment region if the ‘“‘wet adiabatic’’ surface temperature 
(see reference 6 curve in Fig. 12) is less than 32°F. for 
the air speed in question. 

Plotted for comparison purposes in Figs. 9 and 10 are 
the dotted curves of ¢,, for (a) an ice-free surface in clear 
air and (b) an ice-covered surface in clear air. 

In the lower part of Fig. 10 is a plot of the same vari- 
ables at the same altitude but for a lower water-catch 
rate (b = 0.2). Note that the speed at which m = 0 
(ice-free) is almost as high (570 knots) as the corre- 
sponding speed for 6 = 0.5 (587 knots). 

Note by comparing Figs. 9 and 10 that, for a given 
value of } (in this case 0.2), the effect of changing alti- 
tude from 10,000 to 20,000 ft. (without changing the 
ambient temperature) is to lower ¢,, slightly at all speeds 
(except in the 32°F. range), which fact is undoubtedly 
due to the increased rate of evaporation or sublimation 
due in turn to the lower barometric pressure. 


Certain general conclusions can be drawn from a 
study of the four comparative plots. These are: 


(a) High rates of water catch (and correspondingly 
high values of 6) tend to produce higher surface tem- 
peratures in the range of speeds up to about 500 knots. 

(b) In the speed range up to that for m = 1 (dry sur- 
face of ice) for a given value of ambient temperature, 
speed, and 0, the surface temperature decreases slightly 
with altitude. 

(c) The speed beyond which the surface is ice-free, 
even in moderate icing conditions, is much higher than 
the speed for a 32°F. clean surface temperature in clear 
air. For example, in Fig. 10 the speed for » = 0 and 
b = 0.5 is 587 knots compared with 393 knots for a 
32°F. surface at b = 0. 

‘Fig. 11 was plotted to show how this critical ice-free 
speed at = 0 and ¢,, = 32°F. varies with 5, ¢,, and 
altitude. 


(5.3) ‘Wet Adiabatic’’ or ‘Datum Temperature”’ 


As indicated in Part (1), the “datum temperature” 
developed in reference 4 and later used in references 
6, 7, 8, and 9 has been employed in the past as a measure 
of the temperature of an unheated icing surface. The 
following form of this parameter has been copied from 


reference 8: 
V? U? 
ie ces aay: _ Ye | 
to = to + mak a Pr | 


0.622 =(* ) 
Ch P; 
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where 
fo, = datum temperature 
fo = ambient static temperature 
V = free-stream velocity 
U = local velocity just outside the boundary layer 
Pr = Prandtl Number 
L, = latent heat of vaporization (L,) 
€, = saturation vapor pressure of water at fo, 
é; = saturation vapor pressure of water just out- 
side of the boundary layer 
P, = barometric pressure just outside the boundary 
layer 


The term within the bracket represents the local recov- 
ery factor. 

An examination of the above equation indicates that 
this parameter might be more aptly named the ‘‘dy- 
namic wet bulb temperature’ rather than ‘“‘datum tem- 
perature.” Since it does not include terms to account 
for the sensible heating of the impinging droplets or 
the release of the latent heat of fusion, it cannot repre- 
sent the true equilibrium temperature of an unheated 
icing surface or even the temperature of a surface that 
is wetted by an above-freezing-temperature cloud. It 
does represent the temperature in clear air of a surface 
that is wetted by water (possibly from within the sur- 
face) supplied at the temperature of the surface itself. 
This corresponds exactly to the wet bulb temperature 
and, as defined by the above equation, is in a general 
form that can be evaluated at any velocity. In fact, 
it is this form that permits the determination of the 
ambient wet bulb temperature from data obtained in 
high-speed flight. 


Of related interest is the fact that in references 6, 8, 
and 9 the basic equation for the energy transfer from a 
heated airfoil contains two terms that are based on this 
“dynamic wet bulb” temperature. The convection 
loss term is presented as 


Je = felts — tog) 


which would indicate that the difference between sur- 
face temperature and wet bulb temperature is the con- 
This is obviously contrary to the 
If this were actually 


trolling potential. 
definition of the conductance. 
the case, there could be no wet bulb thermometer de- 
pression at low velocity, since the convection gain that 
just balances the evaporation loss would be zero be- 
cause the wick surface would itself be at the wet bulb 
temperature and, hence, the above definition would 
not permit the transfer of heat by convection to the 
wick, 

Similarly, the evaporation term is stated in these 
references as 


de = f(X — 1) (te — tog) 


or 


SCIENCES 


JANUARY, 1958 





0.622L,/e, — ey 
.* 


Cp P, 


Here again, there could be no mass transfer from the 
wet bulb thermometer, since in this case e, would be 
equal to e, and the evaporation term would be zero, 

Fortunately, the two terms q, and g, are invariably 
used in combination with each other in solving anti- 
icing design problems, and, because of the definition of 


the term X, 
a. en - a) 
Cols is _ lo, 


the sum of the above two terms is equivalent to the 
sum of the two corresponding conventional terms: 


YT « 


V2 
c e=Jil, —t hag 
Ye + rd s+ ren) + 


f, (O22) (‘ 7) 
Cp P, 


which satisfy the requirements of the psychrometric 
energy balance. The use individually of either of the 
terms of gq, or q, as defined in references 6, 7, 8, and 9 
would yield inaccurate results. 

In order to illustrate the nature of the divergence 
betweén the results obtained by the methods of the 
present paper and those of references 6, 7, 8, and 9, 
Fig. 12 has been prepared. Curves of ‘datum tem- 
perature” for two ambient temperatures, 0° and 20°F., 
were copied directly from Fig. 2 of reference 6. The 
curves of ¢,, are plotted for an arbitrary value of 0. 
Note that the region of greatest divergence is in the 
speed range below 300 to 400 knots, where the effect 
of the latent heat of fusion is noticeable. In the higher 
speed range, the difference between the two curves is 
small and is apparently due entirely to the cooling ef- 
fect of the impinging liquid which is not taken into ac- 
count in the datum temperature. 


(6) Mopes oF ENERGY TRANSFER TO AN UNHEATED 
IcE-COVERED SURFACE IN CLEAR AIR 


(6.1) t,, Less Than, or Just Equal to, 32°F. 


(6.1.1) Heat Jost by convection: 
Qe = fAlts. — t.) 
(6.1.2) Heat Jost by sublimation: 
gs = 2.90L,f,A[(P;, — P,)/B] 
(6.1.8) Heat gained due to viscous heating: 


Qo = fA(rV,?/2gJc,) 


(6.2) te, Greater Than 32°F. (Wet Ice Surface) 


When the air speed of the ice-covered surface is suf- 
ficient to result in the surface temperature /,, exceeding 
the melting point, a rather complex condition exists. 
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A water film of some indeterminate thickness will cover 
the surface of the ice. The thickness of this film will 
be related to factors such as the liquid-surface tension 
and viscosity, as well as the boundary-layer drag 
forces. The net heat gained due to the viscous heating 
can only produce melting in the ice after having trav- 
ersed the liquid film by conduction. A film of finite 
thickness will therefore contain a finite temperature 
gradient, the differential being (¢,, — 32°F.) since the 
melting ice temperature can only be 32°F. 

The water that is evolved in the melting process may 
leave the surface by two means: (a) evaporation or (b) 
mechanical blowoff or runoff. It has been determined 
by the numerical evaluation of typical operating condi- 
tions that the ice will be melted much more rapidly 
than it can be evaporated, and therefore most of it will 
leave the airfoil by mechanical runoff. For example, in 
order for the rate of evaporation to equal the rate of 
melting at a speed of 800 knots at an altitude of 10,000 
ft. and a conductance of 20 B.t.u./hour °F. ft.?, the 
equilibrium thickness of the water film would have to 
be 0.45 in., an unlikely condition. If it were possible 
for a water film of this thickness to persist, it would 
involve a temperature differential of 31°F. between the 
ice face and the external surface of the water film. If, 
as is probably the case, the actual film thickness is 
small—say, for example, 0.001 in.—its thermal resist- 
ance is small, and it has been determined that little error 
is involved if it is entirely neglected as was the case in 
constructing the curves of Fig. 15. 

The energy balance, based on this latter assumption 
for a wet melting ice surface just slightly above 32°F. 
will therefore include the following terms: 


Heat Jost by convection: 
de = SA (tee a 


.2.2) Heat lost by evaporation: 


(6.2.1) 
t.) = fA(32 — t,) 


— 
ir) 


de = 2.90f-L-A[(Pm — P.)/B] & 


2.90f-L.A [(0.1803 — P,)/B] 


(6.2.38) Heat lost by melting of ice: 
Im = RiumALy 
(6.2.4) Heat gained by viscous heating: 


do = fA (rV.? /2g Icy) 


(7) SUBLIMATION AND MELTING RATES FROM AN 


UNHEATED ICE-COVERED SURFACE IN CLEAR AIR 
(7.1) Sublimation 
The combined heat balance equation for the terms 
listed in Part (6.1) can be written as 
do — Ie = Qs 


or 


UNHEATED 
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V? _— 
felt. + (7p) — | = eaten (F572) 


In the latter form, this equation can be evaluated for 
arbitrary values of ¢,, P,, and B in terms of V, and 
t,. In order to avoid a trial-and-error solution, ¢,, is 
selected as the independent variable. 

Having obtained ¢,, in terms of V,, it is then possible 
to determine the rate of sublimation R,, or y,6,/12 
from 


Qs = 2.9L.f-A[(Ps — P.)/B] = vi(5,/12)L,A 


Since P;, has already been determined, it is possible 
to determine y,6, for an arbitrary series of values of f,. 

Fig. 13 is a plot of such values for three altitudes. 
The corresponding ambient temperatures for each al- 
titude are N.A.C.A. standard. 

The right-hand terminal points of each series of 
curves in Fig. 13 occur at the speed for a 32°F. surface 
temperature, and a general plot corresponding to these 
points is given in Fig. 14. 

Note that in Fig. 13 the rate of sublimation is gener- 
ally small, so that, once ice has accumulated on an air- 
plane, its removal by sublimation alone is a slow proc- 
ess. Note also that for a given value of f, the rate of 
sublimation at 32°F. (right end of each curve) is greater 
at high altitude than at low altitude, principally be- 
cause of the lower value of B and also of P,. Fora 
given speed the removal rate at high altitude is slower 
due entirely to the lower ambient temperature assumed. 
It is significant that these high-altitude removal rates 
would probably require more time for shedding even a 
minor ice accretion than the entire duration of an in- 
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Fic. 15. Rate of melting of ice, 5,4 versus V. in clear air. 


terceptor type of aircraft, the high-speed performance 
of which might be seriously crippled. 


(7.2) Melting 


When the speed exceeds the values given by Fig. 14, 
a wet surface will exist and the controlling factor for ice 
removal is the melting rate, although as indicated in 
Part (6.2) evaporation also occurs but its contribution 
to the removal is small compared to the runoff rate re- 
sulting from the melting. 

The combined heat balance for the terms of para- 
graph (6.2) is 


JANUARY, 1958 


q@-a—- & = de 


or, since /,, = 32°F., 


rV.2 0.1803 — P 
Alt. — = }— 39] —2.9L, 3 
fi E ety | fal B ) 


1 i( bm ‘12)L,A 


This equation was evaluated for the same arbitrary 
series of conditions used in Part (7.1), and the results 
are plotted in Fig. 15. Note that, although the melting 
rates are high, the speeds required to achieve them are 
also high. 

In using Fig. 13 or 15 it should be remembered that, 
as the speed increases, the value of f, also increases, so 
that for a given airfoil the sublimation or melting rate 
will increase faster with speed than indicated by the 
slope of the lines of constant f;,. 


(8) CONCLUSIONS 


At least some ice can collect on an unheated surface 
in normal low-temperature icing conditions at speeds 
up to about 600 knots (690 m.p.h.). 

High rates of water catch tend to maintain an un- 
heated icing surface at 32°F. over a wide range of 
speeds, but, except at the high speed end, a 32°F. sur- 
face temperature does not signify an ice-free surface. 

Surface temperatures predicted by the ‘datum tem- 
perature” or “wet adiabatic’ method of analysis are 
lower in the low-speed range and slightly higher in the 
high-speed range than are obtained by the methods de- 
veloped in this paper. 

Once an unheated surface has been allowed to collect 
ice, the period required to remove it by sublimation 
alone is long. It is of the order of 5 hours for 1/4 in. 
of ice at an air speed of 500 knots at 30,000 ft. and 
N.A.C.A. standard ambient temperature. 

The air speed required to attain a 32°F. ice surface 
temperature in clear air at 30,000 ft. and at N.A.C.A. 
standard temperature is about 850 knots. Beyond 
this speed, melting and evaporation occur. 

There are few reliable experimental data available 
which would permit checking the accuracy of the re- 
sults developed analytically in this paper. Such data, 
obtained in natural icing conditions, would be valuable 
for establishing the validity of the various foregoing 
assumptions (see Appendix). 
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TABLE 1 
EQUILIBRIUM STAGNATION POINT TEMPERATURE OF AN UNHEATED 


LUCITE CYLINDER IN MT. WASHINGTON ICING TUNNEL 


Cylinder Diameter = 
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tunnel static altitude of 7800 feet, 


It was subsequently found that the actual 


average tunnel altitude was more nearly 8200 ft, but the 0 values have not been revised since the net effect on tse 


is believed to be small, 
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APPENDIX 


Since the time of the original presentation of this 
paper, the Lockheed Aircraft Corporation conducted a 
Series of icing tests at the Aeronautical Ice Research 
Laboratory, Mt. Washington, N.H., which facility is 
operated jointly by the U.S. Air Force and the Navy 


Bureau of Aeronautics. One phase of this test program 
included a series of surface temperature measurements 
on an unheated Lucite cylinder, which had a 3.75-in. 
outside diameter and a !/s-in. wall thickness and was 
approximately 12 in. long. Surface thermocouple 
readings were obtained at the forward stagnation 
point, as well as at the 100° and 180° positions. 

The principal purpose of recording these cylinder 
temperatures was to obtain experimental data with 
which to check the validity of the theoretical analysis 
presented in the subject paper. 

The test facility used for this investigation consisted 
of a duct attached to the inlet of a radial-flow jet en- 
gine. This duct, having a cross-sectional area of ap- 
proximately 3 sq.ft., acted as a small icing wind tunnel 
by induction of the ambient icing cloud conditions that 
prevail at the top of Mt. Washington during the winter 
months of the year. By this means, the natural icing 
conditions which flow over the mountain top at veloci- 
ties ranging from about 30 to 90 m.p.h. are accelerated 
in this tunnel to about 250 m.p.h. 

Table 1 contains the data obtained for 18 individual 
icing runs under a variety of meteorological conditions. 
Rotating multiple cylinders were used to obtain the 
liquid water content and mean droplet diameter in the 
vicinity of the test cylinder. The ambient temperature 
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of the air prior to its entry into the tunnel was taken 
to be the most reliable measure of the actual tunnel 
total temperature. The equivalent flight ambient 
temperature was then computed from this total tem- 
perature by making an allowance for the adiabatic ac- 
celeration of the water-laden air, including the effect of 
condensation of additional moisture due to the cooling 
of the air-water mixture resulting from the acceleration. 

From the mean droplet diameter data shown, the 
corresponding droplet Reynolds Number R, and Scale 
modulus y were computed. These values permitted 
the determination of the stagnation point water catch 
ratio (or catch efficiency) Bo and, finally, the theoretical 
value of b, the relative heat factor. The @ values 
shown in the table were determined by using the @ 
curves of the subject paper, together with a set of un- 
published curves which permitted correcting the 0 
values to the altitude indicated in the note under the 
table. 

From the @ values, the theoretical value of n, the 
freezing fraction, and ¢,,, the equilibrium stagnation 
point temperature, were determined. In the last 
column the actual measured values are tabulated for 
comparison. The agreement is extremely good in spite 
of the obvious limitations of the instrumentation. For 
example, in Run No. 2A, a measured temperature of 
34°F. is indicated even though the cylinder is unheated, 
an ice formation was observed in the thermocouple 
region, and the tofa/ temperature of the tunnel was 13.2° 
F. Similarly, the measured temperatures of 33°F. in 
runs 3, 8, and 9 indicate that the thermocouple inac- 


curacy was as high as 1 or 2°F. under some conditions, 
It is interesting to note that, in runs 5 and 13 where the 
predicted values of ¢,, are 1° or 2° higher than the meas- 
ured, the value of 1 is only slightly less than 1.0, which 
would be the true value of m based on the measured 
temperature. 

Note in the second column that the proportion of the 
runs obtained in natural icing was about equal to that 
using artificial spray water icing. The fact that the 
agreement in the last two columns did not appear to be 
influenced by this factor is an indication that, under 
the proper control, artificial icing is a valid research 
technique. The fact that it has been reported to the 
writer that the full effect of the latent heat of fusion has 
not been observed in the N.A.C.A. icing tunnel or in 
the Canadian N.R.C. icing tunnel is an indication that 
the spray water systems used in these facilities are not 
producing 100 per cent subcooled liquid droplets and 
may be hampered by the presence of partial “‘freeze- 
out” prior to impingement. That this does not seem 
to occur in the Mt. Washington tunnel may be the for- 
tuitous result of the necessarily short distance, and, 
hence, time span, between the spray nozzles and the 
test section. 

In conclusion, it is the opinion of the author that the 
experimental data presented in this appendix are strong 
evidence that the basic equations of the foregoing paper 
are valid. The extent of the agreement shown in the 
table between measured and predicted temperatures 
cannot be mere coincidence, since the conditions of the 
tests covered a wide range of icing conditions. 
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Blade Adjustment in Axial-Flow 
Compressor Stages’ 


ALLAN W. McCOY? ann ALTON V. HOOPER? 


Sverdrup e& Parcel, Inc.** 


SUMMARY 


An investigation of blade adjustment as a means of achieving 
volume regulation of an axial-flow compressor stage is presented. 
The effectiveness of such blade control is established and defined 
as the relative change of volume flow per unit of adjustment angle 
at constant pressure rise. 

In the case of a grid element, the effectiveness of blade control 
is shown to be dependent upon the principal design parameters of 
reaction ratio, flow, and pressure coefficients. Low reaction 
ratios should be employed for stator-blade control, and for corre- 
sponding use of rotor-blade adjustment the preference is toward 
high reaction ratios. 

Radial distribution of blade-control effectiveness is shown for 
several conventional stage types without consideration of radial 
equilibrium and nonuniform energy addition. In case of stator 
grid adjustment, the distribution of effectiveness indicates maxi- 
mum values at the blade root section—i.e., with the adjustment 
for increase in flow, the axial velocities tend to rise rapidly at the 
hub. Rotor-blade control in these conventional stage diagrams 
affords a somewhat more uniform distribution of effectiveness. 
Special stage designs for uniform spanwise distribution of adjust- 


ment effectiveness are also discussed. 


INTRODUCTION 


F )R EFFICIENT COMPRESSION OF FLUID with high 
rates of flow, the axial-flow compressor offers 
definite advantages over competitive turbomachines. 
A significant limitation, however, is its narrow range of 
flow at constant speed due to the inherent steepness of 
the pressure-volume characteristic. Certain applica- 
tions of axial-flow compressors require a wide range of 
both pressure ratio and volume flow, and the usual 
compressor characteristics must be modified in some 
manner. ‘This requirement of a wide operational range 
is particularly evident in the case of compressors for 
high-speed wind tunnels, such as those at the Arnold 
Engineering Development Center of the U.S. Air Force. 
In such high-speed wind tunnels, the pressure and flow 
requirements are no longer satisfied by the usual com- 
pressor performance characteristics with variable drive 
speed; anew method of effective wind tunnel and com- 
pressor matching must be employed. Because of the 
large size of these wind tunnels and compressors, con- 
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siderable benefits accrue from application of design 
principles permitting high mass-flow ratios! and attain- 
ment of required operating range by use of adjustable 
blading at a constant rotational speed. 

Although the use of stator-blade control offers less 
mechanical difficulty than similar adjustment of the 
rotor blades, both types of control will be considered for 
purposes of generality and comparison. Considera- 
tions are based upon a single “normal’’ stage that is 
capable of multistage application. 


NOTATION 


A listing of the nomenclature is given below. In addition, a 
sketch of a typical blade diagram of a normal stage, with much 
of this notation indicated, is shown in Fig. 1. This normal stage 
velocity diagram is such as to produce exit flow conditions corre- 
sponding to those prevailing at the stage inlet. The major por- 
tion of the notation used in the subsequent analysis is given by 
Rannie.? The reference velocity for securing this nondimension- 
ality is usually the peripheral velocity of the rotor tip, although in 
some cases the rotational velocity of a particular element is used. 
In this latter instance, all flow terms are said to be ‘“‘local’’ values. 


= flow coefficient (dimensionless axial velocity) 


rr) 

A = rotational component of absolute velocity 

“ = rotational component of relative velocity 

y = pressure coefficient (dimensionless pressure rise) 

U rotational speed of rotor tip element, used for nondi- 
mensionalizing above velocity terms 

~ = radius ratio 

R = reaction ratio (ratio of static pressure rise in rotor to 
static pressure rise of complete stage) 

n = relative change of compressor flow or axial velocity 

E = effectiveness of blade control 

8 = rotor flow angle 

y = stator flow angle 

6 = reset angle of adjusted grid 














nll 


= . 
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Fic. 1. Typical compressor diagram with stator grid reset. 
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Subscripts 
R = rotor 
S = stator 
0 = hub conditions 
1 = rotor inlet conditions 
2 = rotor exit conditions 


GENERAL THEORY 


(1) Definition of Control Effectiveness 


The use of blade adjustment in axial-flow compressor 
stages results in relocation of the constant speed pres- 
sure-volume characteristics. Such typical characteris- 
tics are indicated in Fig. 2 for resetting of either the 
stator or rotor grid. In this case, as well as in the fol- 
lowing sections, the grid adjustment angle is taken as 
positive when the action is such as to increase the vol- 
ume flow. 

For establishing the effectiveness of blade control, 
attention is centered upon design operation at zero 


resetting and operation at an increased flow with ad-. 


justment of either rotor or stator grid by an amount 
6. When either grid is reset to this latter position, an 
operating condition 


te = o+ Adel 5 =a (1) 
Ve=vt AVS xv 


exists where the blade grids are aerodynamically loaded 
equivalent to design values. In general, this equiva- 
lent operating point would be at a slightly different 
pressure rise than the original value. However, blade 
adjustment is considered primarily as a means of vol- 
ume regulation, since the usual multistage compressor 
characteristic approaches constant volume operation. 
The new operating point of particular interest is at 


6 = + Adl , _ 
v=" J 


(2) 





Sb 











> 


Typical constant-speed characteristics with resetting 
of rotor or stator grid. 


Fic. 2. 


1953 


This latter condition, representing a change of flow at 
constant pressure rise, provides a criterion for defining 
the blade-control effectiveness. 

Let the relative increase of flow be given by 


n = A¢/o = (¢’ — $)/¢ (3) 


at y’ = y. Then the rate of change of relative flow 
with respect to adjustment angle under these conditions 
is 


(On/05), = n/éb (4) 


which will depend upon the blade diagram design and 
the amount of grid adjustment angle. In order to 
compare different types of blading design, this deriva- 
tive may be evaluated at the original design conditions 

i.e.,at 6 = 0. It then represents the initial slope of 
the n — 6 tunction at constant stage pressure rise and 
provides a definition of the blade-control effectiveness. 
Thus, 


E = (On, 05) y, até =0 (5) 


(2) Assumptions 


For any analysis of compressor blading, certain con- 
ditions must be imposed. The assumptions that are 
presently introduced are: (1) incompressible fluid; 
(2) radial flow components in the increased flow ob- 
tained by blade adjustment are similar to those exist- 
ing at design conditions for each blade element; (3) 
constant fluid exit angles from grid elements; and (4) 
uniform radial distribution of stage pressure rise. 

The assumption of incompressible fluid, while perhaps 
compressor,’’ is common to 


“cc 


inconsistent with the term 
most blade theories for purposes of at least first approxi- 
mations. Methods of compressibility correction for 
higher order approximations exist. 

The second of the above assumptions is not intended 
to exclude any possibilities of any radial flow. It is 
intended to admit the existence of radial flow compo- 
nents at the design condition, determinable by radial 
equilibrium conditions,*~* but assumes that similar 
radial flows exist in the off-design flow obtained by blade 
adjustments. It is expected that this condition may be 
rather closely realized when the grid reset angles are 
small, and in this case the effectiveness of blade ad- 
justment is evaluated at the design conditions with 
zero reset. 

The assumption of constant fluid exit angles is neces- 
sarily introduced in view of the present inability to 
specify variations of deviation angles except in terms of 
a particular grid configuration. For more exact calcu- 
lations, the known characteristics of the particular cas- 
cade considered may be substituted for this assumption. 

While the assumption of uniform rotor work in a 
normal stage is common to most direct theories, it is 
strictly valid only at the stage design conditions. For 
off-design operating conditions, this assumption consti- 
tutes somewhat of an oversimplification and leads to 
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only first-order determination of the flow distributions 
at the larger departures from design flow. However, 
uniform energy addition is physically desirable for 
multistaging purposes. The flow distributions ob- 
tained by this assumption may be employed as the basis 
of further approximation by direct numerical methods 
with variable energy rise. The magnitude of induced 
effects and radial variation of pressure rise may be gen- 
erally estimated as proportional to nonuniformity of 
blade-adjustment effectiveness, as determined by the 


present method. 


(3) Stator Blade Adjustment 

For a fixed position of the rotor grid, the change in 
volume flow may be determined in terms of the stator- 
blade adjustment angle, ds. 

At a constant pressure rise, Y = y’ and 


AX = Ad’ = y/2é (6) 


where the primes denote values obtained at stator reset 
angle, ds. Radial equilibrium effects may cause the 


condition that 
oi ~ od» 


at a particular radius. 

Therefore, let the change of volume flow rate (axial 
velocity) at an arbitrary radius be given in terms of both 
the rotor inlet and exit conditions, 


oi’ — 1 go’ — 2 
od | : o: 


mm = 


The condition that the relative flow exit angle remain 
constant may be written for the rotor as 


(E — do’)/p2’ = (E — Ao)/de (8) 


and for the stator grid, 
(o:’/M1’) — (b1/Mi) 


pronase ; (9) 
1 + (1/1) (G1 /As) 


tan 6 = tan (y1' — vi) = 


Eq. (8) may be written in terms of the rotor inlet 


conditions as 


2’ 


— ri’ — (9/28) 


: " 
| &—m — (¥/28) (Sa) 


From the assumption that the increased flow produces 
similar radial velocity components at each radial ele- 
ment, 2m; = mo = n, and 
g2'/b2 = gi'/b1 (10) 
Introducing Eqs. (Sa) and (10) into Eq. (9), the sta- 
tor-adjustment angle is given by 
p 
a (9a) 
nN a. A? + gi- 
n §& — (/2&) 


The relative flow variations may be obtained as 


tan 6 = 
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Fic. 3. Blade control effectiveness versus reaction ratio for 
grid element. 
n = B/(A cot6+ 1 — B) (11) 
where 
a-% pot G@rt/r) 
ny — — (¥/28) 


Accordingly, the derivative of the m — 6 functional rela- 
tion 1s 
(5) _" B sec? bs 4 
05s/ y A+ (1 — B) tan ds 
B tan 6s(1 — B) sec? dbs 


12 
[A + (1 — B) tan ds]? (12) 


and the stator-blade-control effectiveness, by definition 


established in Eq. (5), is 
A? + 1” 


Ee = (m) — : (13) 
“ N\Ws/ys-0 A dle — (¥/28)) °° 


(4) Rotor-Blade Adjustment 

The effectiveness of rotor-blade control with a fixed 
position of the stator grid may be similarly determined. 
In this case, the constancy of the stator grid exit flow 


angle is 
A’/ or’ = Ai/hr (14) 


and the corresponding condition for the rotor blades 
with adjustment angle, 6g, may be expressed as 


, / / 
tan dr = tan (B2.’ — B2) = (2/2) — (oe He) (15) 
1 + (o2’/p2") (b2/ me) 
where pe = & — dz. 
By a procedure similar to that used in the case of 
stator-blade adjustment, the effectiveness of rotor- 
blade control may be obtained in the form 


On M2” + do»? 


Ep = ) = | 16) 
. fs vde=0 G2l& — (¥/28)] ane 
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EFFECTIVENESS OF ELEMENT 


To permit comparison of rotor and stator-blade con- 
trol, it is desirable to consider the flow conditions at a 
particular radius. In this case, an arbitrary radius 
ratio, ~*, was selected in such a manner that 


d=¢ = (até = &) (17) 


With radial equilibrium this condition is usually ob- 
tained at essentially mean radius, and the degree of re- 
action is defined simply as 


R = 1 — (An/€*) (18) 


where Amn = Ai + (AX/2). 

Since attention is now confined to this particular 
grid element, the reference velocity may be changed 
from the rotational speed at the rotor tip to that at the 
particular element. Thus, in terms of ‘“‘local’’ values, 


a {1 — R — (¥/4) + ¢ 
o[1 — (¥/2)] 
IR — W/4)P + # 

o[1 — (¥/2)] 
By these equations, the effectiveness of both rotor and 
stator adjustment of a grid element is expressed in terms 
of the fundamental design parameters of reaction ratio 
and pressure and flow coefficients. The independent 
influence of these design parameters is shown in Figs. 
3-5. 


Es 





(19) 








Er = (20) 


Reaction ratio effects are shown in Fig. 3 for particu- 
lar values of flow and pressure coefficients. From Eqs. 
(19) and (20), it may be noted that the reaction ratio, 
R, is the only design parameter that gives rise to differ- 
ences between rotor- and stator-blade-control effective- 
ness. Furthermore, both types of adjustment have 
equal effectiveness at 50 per cent reaction ratio. As 














Blade control effectiveness versus flow coefficient for grid element. 


indicated in Fig. 3, with increase of reaction ratio, 
the rotor-blade-control 
whereas the corresponding effectiveness of  stator- 
blade-control decreases. 


effectiveness is increased, 


The effect of flow coefficient on the effectiveness of 
rotor- or stator-blade adjustment is shown in Fig. 4. 
These curves indicate the general desirability of high 
flow and pressure coefficients in attaining high effective- 
ness of blade control. The difference between rotor- 
and stator-blade control is dependent only upon the 
reaction ratio, as noted previously. 

The variation of effectiveness with pressure coef- 
ficient is seen to be dependent upon the values of the 
other design parameters. Certain further insight may 
be had by examining the conditions for the mini- 
mum effectiveness values shown in the E — ¢ plots of 
Fig. 4. Differentiating Eq. (19) for the stator-blade- 
control effectiveness with respect to the flow coefficient, 

OF s 1 A? 





es —— — ——*_____ (a1) 
a = 1-—(¥/2) gl — (¥/2)] 


and the condition for minimum effectiveness of stator- 
blade control is that 

1 =¢ (22) 
which occurs at a rotor inlet flow angle of 7/4. The re- 


sulting expression for the minimum stator-blade-con- 
trol effectiveness is 


E ee. eo 1 — 2R- 93) 
wae" 1- W/D 1-w2 “ 


The corresponding condition for minimum rotor-con- 
trol effectiveness is 


Me = (24) 


and the minimum value is given by 
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Ermin. = 1 — {(1 — 2R)/[1 — (¥/2)]} (25) 


Again, these relations for rotor and stator control are 
symmetrical about 50 per cent reaction and are shown 
in Fig. 5. It is indicated that the advantages of in- 
creased pressure coefficient depend upon which control 
method is employed and whether the reaction ratio is 
greater or less than one-half. The minimum effective- 
ness of stator-blade control increases with pressure coef- 
ficient provided the reaction ratio is less than one-half. 
This condition for rotor-blade control requires reactions 
of greater than one-half. 


BLADE CONTROL IN COMPLETE STAGES 


The relations for blade-control effectiveness have 
been given in Eqs. (13) and (16) for stator- and rotor- 
blade adjustment and were expressed in terms of the 
relative rotor inlet and absolute rotor discharge condi- 
tions, respectively. By use of a simplifying assump- 
tion that 


oi = 2 = > (26) 


at all radii, it is possible to make immediate compari- 
sons of various stage designs for application to blade 
control. For this comparison the equations for blade- 
control effectiveness are 


Es = (A? + ¢°)/o[& — (/28)] (13a) 
Er = (u2> + ¢°)/o[— — (¥/28)] (16a) 


Several conventional compressor-stage types may be 
specified by writing the absolute compressor inlet spin 
velocity as 


\ = Ke* (27) 


and the flow types obtained with X = —1, 0, and +1 
will be considered. The blade diagrams corresponding 
to these selections are of fundamental nature and are 
frequently discussed for other comparative purposes. 
By conventional reference they shall be called: 

Free Vortex (FV), X = —1 

Constant Inlet Spin (CIS), X = 0 

Wheel Flow (WF), X = 1 


The constant K may be selected according to the reac- 
tion ratio desired at a particular radius. 


It is also considered desirable to include comparison 
of blading types that would offer uniform effectiveness 
distribution along the blade span. By definition, these 
uniform effectiveness distribution blading types would 
satisfy the relation 


OE /0E = 0°n/O0EOS = 0 (28) 

Description of the flow in terms of the absolute rotor 

inlet spin velocity may be made for comparison with 
conventional blade types given by Eq. (27). 

The velocity diagram for uniform distribution of sta- 

tor blade adjustment effectiveness (USED) is given by 
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d= {Esolt — (W/2E)] — 6?}'7 (29) the conventional-stage types considered. In some 


and for uniform rotor control effectiveness distribution 
(URED) by 


1/2 
MW = E- : - E (« ~ ) ~ é| (30) 


where Eg and Fp are selected constants. 

Plots of the spanwise effectiveness distributions for 
both rotor- and stator-blade control with the stage types 
considered are shown in Figs. 6 and 7. The rotor inlet 
spin components are also indicated for comparisons. 
In Fig. 6, the appropriate constants have been selected 
such as to produce a symmetrical velocity diagram at 
the blade root section. The overall effectiveness of 
these stages is seen to vary appreciably, and a further 
comparison, with a symmetrical diagram at mean ra- 
dius, is shown in Fig. 7. 

The examples of blade control in conventional axial- 
flow compressor stages show that the control effective- 
ness is, in general, nonuniformly distributed over the 
blade span. For wide departure from design condi- 
tions, the flow pattern will be influenced by such condi- 
tions as radial unbalance and possible secondary flow 
development. These effects will be most predominant 
in blade types that exhibit highly nonuniform distribu- 
tion of blade control effectiveness here. Certain addi- 
tional qualitative information may be deduced from 
the investigations by application of simple radial un- 
balance considerations. 

All cases of stator-blade control, shown in Fig. 7 for 
conventional diagrams, indicate an effectiveness dis- 
tribution such that maximum values are obtained at the 
blade root sections. Thus, for an increase of flow, the 
axial velocity profiles will tend to peak at the hub and 
reduce the aerodynamic loading of these root elements. 
Consideration of radial equilibrium effects indicates 
that the flow will then readjust itself such that these 
excess axial velocities will be reduced somewhat and 
the stage energy addition will contain some positive 
radial gradient. This positive energy gradient is in 
the direction of increased flow stability.° Similar 
stator-blade adjustment for decreased flow will produce 
an opposite effect, leading to flow starvation at the hub 
elements, negative radial energy gradient, and possible 
instability. These effects may persist and result in 
reduction of the attainable pressure rise and operating 
range at less than design flow. 

Rotor-blade adjustment is seen to produce a more 
nearly uniform distribution of blade-control effective- 
ness than corresponding stator-blade adjustment for 


cases, stator and rotor control complement one another 
to some extent, and a method of combination contro} 
could be employed to provide an essentially uniform 
distribution of the combined blade-control effectiveness, 

While both types of blading for uniform blade-control 
effectiveness differ from those in conventional use, the 
required flow patterns do not appear unduly radical, 
This is particularly true of the uniform rotor-control 
diagram, which appears to resemble a combination of 
free-vortex and wheel flow characteristics. Since radial 
unbalance and limit loading conditions were not con- 
sidered, further comparison of stage types should in- 
clude such effects. 


CONCLUSIONS 


From considerations of the blade-control effective- 
ness of a grid element, it is indicated that high effec- 
tiveness values may be attained by use of high values 
of flow and pressure coefficients. Volume control by 
either rotor- or stator-blade adjustment produces sym- 
metrical behavior of the control effectiveness about a 
50 per cent reaction ratio axis. Low reaction ratios 
should be used for stator-blade control and high values 
for corresponding rotor adjustment. 

The radial distribution of blade-control effectiveness 
in several conventional-stage diagrams is shown to be 
generally nonuniform. In all cases of stator-blade 
control and in some instances with rotor-blade control, 
the effectiveness is a maximum at the blade roots. 
This spanwise distribution is somewhat more uniform 
in the case of rotor-blade control. Special stage types 
for uniform distribution of stator- or rotor-blade control 
effectiveness do not appear unduly extreme and further 
analysis to include radial unbalance and loading limita- 
tions is suggested. 
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Generalized Theory of Convective Heat 


Transfer in a Free-Molecule Flow’ 


A. K. OPPENHEIM?t 
Unwersity of California at Berkeley 


ABSTRACT 


The theoretical analysis of convective heat transfer in a free- 
molecule flow with Maxwellian velocity distribution is generalized 
by reducing all relationships to a form common to all gases irre- 
spective of their molecular structure. A simplified c®mputing 
procedure is set up, and the heat-transfer characteristics of a flat 


plate, cylinder, and sphere are presented and correlated graphi- 


cally. 


A method of applying these results to a cone and to com- 


pound shapes is described. 


Io(s) 


Ii(s) 


lerfcs 


Nu 


NOMENCLATURE 


specific heat capacity at constant pressure, 
(ft.lbs.)/(Ib.sec.2/ft.)°F, 

specific heat capacity at constant volume, 
(ft.lbs.)/(1b.sec.?/ft.)°F. 

reference diameter, (ft.) 

molecular energy transport rate, (ft.lbs. /sec. ) 

molecular energy transport rate per unit surface area, 
(ft.lbs. /ft.?sec. ) 

base of natural logarithms 2.7183 


9 s 

error function, e~2? dx, (dimensionless) 
VF Jo 

complementary error functions, 1 — erfs, (dimension- 


less ) 


average surface integral defined by Eq. (2.13), 





f dS, (see Q) 
ws 
average surface integral defined by Eq. (2.14), 


g dS, (see Q) 
S 
heat-transfer coefficient, Q/(Tw — Taw), 
(ft.lbs./ft.2sec. °F.) 
modified Bessel function of first kind and zeroth 


af 


T 
order, e— 5 ©°S ® dg, (dimensionless ) 
0 
modified Bessel function of first kind and first order, 
1 7 


_ e~ 5 ©°S % cos @ dd, (dimensionless ) 
¥ Je 


integral of the complementary error function, 


erfcx dx, (dimensionless ) 
Ss 
number of degrees of freedom, (dimensionless) 
Boltzman constant 5.66.10~24 (ft.lb./molecule°F.) 
Mach Number (dimensionless ) 
number of molecules per unit volume (molecules /ft.*) 
Nusset Number 4D/x (dimensionless) 


Received May 21, 1952. 
*This work was done as part of an Office of Naval Research 


contract with the Department of Engineering, University of 
California at Berkeley. 
} Assistant Professor of Mechanical Engineering. 


ll 


ll 


ll 


ll 


number of molecules striking unit surface area per 
unit time, (molecules/ft.*sec. ) 

Prandtl Number, cpu/x (dimensionless ) 

heat-transfer rate (ft.lbs./sec.) 


; ] 
average heat-transfer rate per unit area, 5 q dS 
' Ss 


(e.g., for a cylinder Q¢ = q 40; for a sphere 


9 
—a/<é 


1 w/2 
Os = q cos 6 d@), (ft. lbs./ft.? sec.) 


“J —x/2 

local heat-transfer rate per unit area, (ft.lbs./ft.*sec.) 

Reynolds Number, D U’p/y, (dimensionless ) 

recovery factor, (Taw — To)/(Ts — To), (dimension- 
less) 

surface area, (sq. ft.) 

Stanton Number, h/cp,U = Nu/RePr, (dimension- 
less ) 

molecular speed ratio, U/ Vm = V/ y/2 M, (dimension- 
less ) 

temperature, °R. 

mass velocity, (ft. per sec.) 

internal energy, (ft lbs./molecule) 

most probable molecular speed, (ft. per sec. 

accommodation coefficient, (e; — er/(e; — ew), (di- 
mensionless ) 

Cp/Cy, (dimensionless ) 

U/Vm = s sin 6, (dimensionless) 

angle of incidence, (deg.) 

thermal conductivity, [ft.lbs./(ft.2sec. °F./ft.)] 

viscosity, (Ib. sec./ft.?) 

density, (Ibs.sec.?/ft.*) 

(Taw/To) — 1, (dimensionless ) 

(NVm 2/x \(e~”"/2) (molecules /ft.? sec.) 

e-” + (3/2)Vxnl(1 + erfn), 
e-™ + y/an(1 + erfn) 


(dimensionless) 


: + 


adiabatic wall (equilibrium) condition 

front side of body with respect to flow direction 
incident molecules 

flat plate normal to flow direction 

free-stream condition 

rear side of body with respect to flow direction 
re-emitted molecules 

stagnation condition 

the translational ene1 gy 

(incident) molecules at wall temperature condition 
velocity component normal to surface 

flat plate at an angle of incidence 


(1) INTRODUCTION 


HE FREE-MOLECULE FLOW REGIME represents one 
the few instances, besides some examples 
in laminar flow, where the convective heat-transfer 
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characteristics can be deduced entirely from funda- 
mental principles. To date, the expression for the heat- 
transfer rate to a plane surface exposed toa free-molecule 
stream with Maxwellian velocity distribution has been 
derived,! and the results integrated for, transverse flow 
past a circular cylinder? and for a sphere.* Further- 
more, a satisfactory experimental check for the case of 
the cylinder has been demonstrated.‘ 

In all the above-mentioned analyses, separate expres- 
sions were derived for a monatomic and a diatomic gas. 
While such a distinction is obviously necessary for 
practical applications, it does not have to be carried 
through in the analysis where the gas properties depend- 
ing on the molecular structure can be represented by a 
suitable function of the number of degrees of freedom, /, 
or of the ratio of specific heats, y. The results then be- 
come general as far as the molecular structure of the 
substance is concerned. Furthermore, by correlating 
the results for a few fundamental surface shapes, a basis 
is formed for a simple method of evaluating the heat- 
transfer characteristics for a great variety of surfaces. 
This is due to the fact that the results are functions of 
average values of integrals taken over the surface for 
which the heat-transfer characteristics are evaluated, and 
the integrands are the same for all surfaces. Con- 
sequently, the computation for any compound surface 
will consist in simply evaluating the average solution 
for the component fundamental shapes. 

In this paper the theory of the convective heat trans- 
fer in a free-molecule flow in Maxwellian equilibrium is 
generalized and correlated in the sense indicated in the 
preceding paragraph. The results are applicable to a 
substance of any molecular structure and cover the 
range of all the extreme and simple surface shapes. 


(2) ANALYSIS 
(2.1) Energy Equation 


The convective heat transfer from a body in a free- 
molecular flow is governed by the energy balance, 


Gq= CC, — & (2.1) 


where g = dQ/dS is the convective heat-transfer rate to 
the gas stream per unit surface area of the body, e; = 
dE,/dS is the energy transport rate of the incident 
molecules per unit surface area, and e, = dE,/dS is the 
energy transport rate of re-emitted molecules per unit 
surface area of the body. 

Since the number of re-emitted molecules depends on 
the “‘accommodating”’ property of the surface, which can 
be determined only by experiment, the energy associ- 
ated with these molecules is usually expressed by means 
of the accommodation coefficient defined as 


a = (e; — e,)/(€; — ew) (2.2) 


where e, =dE,,/dS and E, is the flow rate of the energy 
that would be reflected from the surface if all the mole- 
cules were re-emitted with Maxwellian distribution at 
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the body wall temperature. Eq. (2.1) then becomes 


g/a = Cy — & (2.3) 


Each of the energy terms appearing in the above 
equation consists of the translational energy plus that 
part of the internal energy which is due to other de- 
grees of freedom than those associated with transla- 
tional motion. 

By the principle of equipartition, the internal energy 
of a gas molecule with j degrees of freedom in equi- 
librium at temperature 7'is u = (j/2)kT, where k is the 
Boltzman constant, and, since 


fp _ (/aT)utkT) _ |, 2 


Ce (du/dT) 7 j 


ll 


y 
od 


one obtains 
u = [1/(y7 — 1)JAT. 


Now, the translational energy of a molecule issuing 
from the wall at temperature 7% is @4. = 2nkT,,—that 
is, 4/3 times as great as the mean translational energy 
of the gasin equilibrium at the same temperature.’ Since 
the translational part of the internal energy of a gas in 
equilibrium is due to three degrees of freedom, it follows 
that 


— 3 j 1 
Cw = Cw os > nkT — i a nkT = 
l 
yt) er. (24) 
Ay — 1 


where n is the number of molecules striking unit sur- 
face area of the wall per unit time. 

Similarly, the energy transport rate of the incident 
molecules is 


— 3 
nkTy = e - 
2 ae 1) 


J . , 
e&=e@+ —— nkTy (2.5) 
where e, = dE,/dS is the translational energy flow rate 
of the stream striking a unit surface area of the wall per 
unit time and 7» is the free-stream temperature. 


Introducing Eqs. (2.4) and (2.5) into Eq. (2.3) results 
in 
| i+ 


1 —s 
ie 9 nkT» — («. += nk Ts) (2.6) 


Qa “a 
or, in terms of y, 


1 y+1 ; 5 — 37 r| 
at sh. sape pennienncierints Rl» _ ' _ kT» 
. t-0" jot 


which, with the proper values of 7 = 3 or y = 5/3 fora 
monatomic gas and 7 = 5or y = 7/5 for a diatomic gas, 


7 
reduce readily to Eqs. (30) and (31) of reference | or 
Eqs. (1) and (2) of reference 3. 








(2.2) 


give 
stril 
flow 
into 
surf 


flow 
mol 


shor 


whi 
by\ 
exp 
whi 
fun 


et, 


or, 


ank 


mo: 
ave 


tio 


al 





mes 


(2.3) 


above 
> that 
‘r de- 
insla- 


ergy 
equi- 
is the 


wing 
that 
ergy 
ince 
1S inl 
lows 


2.4) 


sur- 


ent 


ate 
per 


Its 








CONVECTIVE HEAT TRANSFER 


(2.2) Average Heat-Transfer Rate 


In order to determine the heat-transfer rate from a 
given body, expressions for the number of molecules 
striking a unit surface area per unit time, #, and for the 
flow rate of translational energy, e,, must be introduced 
into Eq. (2.6) and the equation integrated over the body 
surface. The expressions are as follows: 


NVn 
2/ a 
e, = {[s? + (5/2)]n — O}RT 


le + Va n(1 + erfn)] (2.8) 


a = 


(2.9) 


where VV is the number of molecules per unit volume; 
Vm is the most probable molecular speed; 7 = U,/ Vm, 
U, being the component normal to wall of the mass 
flow velocity U; s= U/Vn = Vy 2 M is the so-called 
molecular speed ratio, which is related directly, as 
shown, to the Mach Number, 1/7; and 


© = (NV 4o/ me” (2.10) 


The derivation of Eq. (2.8) is described in reference 6, 
while an equation equivalent to Eq. (2.9) was obtained 
by Stalder and Jukoff in reference 1. In the latter, e; is 
expressed in terms of a function y (see nomenclature), 
which has been here replaced by a somewhat simpler 
function ®, defined by Eq. (2.10). 

By the substitution of Eqs. (2.8) and (2.9) for m and 
e;, Eq. (2.7) becomes 


] I 
re = nkT -|( +— - ) n — »| kT 
a 2(y — 1) y¥-—1 


-g= 
2.11) 
or, in a nondimensional form, 


q 77 1 n Tv 


ankT)U = 2(y — 1) NV ns To 


( ‘e— ) foo (19) 
= 7 . ole 
» ices 1 N J se N J ee 


The integration of the above equation is performed 
most conveniently by evaluating first the following 
average surface integrals: 


] e” 
C=- - dS 
Lf 2/ as 


is uff n(1 + erfn) ds 
Ses 2s 


S denoting the surface area of the body under ques- 
tion. Then Eq. (2.12) yields 
Q ¥ +1 a, I, 
= (G+ Fh) - 
Ay — 1) To 


2 


2.13) 


(2.14) 


aNkT,U 
Bi antl 
+) G+F)+<=G (2.15) 


where 


_ 


FREE-MOLECULE FLOW 5 


I 
= fads 


is the average heat-transfer rate by convection per unit 
surface area. It might be of interest here to note that 
Gand F are related to the function x, used in reference 1, 


as follows: 
= l x . 
G+F=-— z= S$ 
SSS 2 xs 


(2.3) Heat Transfer Parameters 


IN A 


OI 


All the essential results of the heat-transfer analysis 
in the free-molecule flow are most conveniently repre- 
sented in terms of the Stanton Number, St, and the re- 
covery factor, r. It is of interest to note that the rela- 
tionships between these parameters and the Mach 
Number, /, can be again represented in a general man- 
ner for any gas independently of its molecular structure. 

For 0 = 0, Eq. (2.15) yields the expression for the 
adiabatic wall (equilibrium) temperature to the free- 
stream temperature, 


‘em mic I I 7 
Tew _ 9% ( ae :) 2.16) * 
T» TT, TVW 1 2G+H 


Introducing the above into Eq. (2.15) results in 


= > — er - 
Gna (G+ F)(T, — Ton)NRU (2.17) 
2(y — 1) 
and the Stanton Number 
h ) + 1 " 
St = — == —— = ee. a(G + F) 
pCyl (7, — Ten) tol 2 
(2.18) 
since 
yt+1 Nk vyt+tl1q-q vt! 
Ay — 1) pe, 7 Ay — 1) Sy 7 2y 
It should also be noted that 
Nu 9y — 5 Nu 
Hn ae Sn (2.19) 


ss RePr 4y_ Re 


since, in the free-molecule régime, Eucken’s expression 
for the Prandtl Number, Pr, in terms of y is well justi- 
fied. With Eq. (2.19), therefore, Eq. (2.18) can also be 
regarded as expressing the relationship between the 
Nusselt Number, Nu, the Reynolds Number, Re, and 
the Mach Number, 1. The latter does not appear ex- 
plicitly in Eq. (2.18), but, it should be remembered, the 
values of G and F are, for a given body geometry, only 
functions of the molecular speed ratio, s, or Mach Num- 
ber. 

Finally, by the use of Eq. (2.16), the recovery factor 
can be expressed as 


* In integral form valid only, of course, if 7» is constant over 


the body surface—i.e., if the body conductivity is infinite 
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since the ratio of the stagnation temperature 7), to the 





free-stream temperature 7) is 
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It is obvious from Eqs. (2.18) and (2.20) that, in 
order to represent the results of the analysis in a general 
form independent of the molecular structure of the sub- 
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. a -03 == SS SS SS SS ES SSS 
stance, one should express the results in terms of 04 tt =a} || | 
Ria ) tq a , oe. ; > . 7 ss -05 oe i it i 
ly/(y¥ + 1) |St instead of St and in terms of [(y + 1) + 68s ea oe ert SS 4 a 
y|rinstead of r. It is of some interest to note here that Mop 03 040806 08 i ¢ a mae Ge wr 
the product of the Stanton Number and the recovery Mz02 03 040506 0810 15 2 OT "tt > & 


Fic. 2. Average surface integral F as a function of the mo- 


factor is independent of the molecular structure of the L s 
lecular speed ratio (and the Mach Number) for a flat plate, a 


substance. sphere, and a transverse circular cylinder. 
(3) DISCUSSION to. 
a St2 = gr7 ot 
(3.1) Results ‘Tol: Sik: Sac Ge | T vonnety anny 
ee | r og] SUBSCRIPT | REFERS TO A MONATOMIC GA | 
4} t SUBSCRIPT 2 REFERS TO A DIATOMIC GAS | 


According to the method developed in the preceding 
section, the analysis of the convective heat transfer in a 
free-molecular flow consists of computing first the values 
of G and F defined by Eqs. (2.13) and (2.14) and then 


determining 





* In integral form valid only, of course, if 7, is constant over = Sata ee 
the body surface—i.e., if the body conductivity is infinite 
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SUBSCRIPT | REFERS TO A MONATOMIC GAS | M2 02 03 04 e a 5 2 » = 6 8 M2 
se — Fic. 3. Stanton Number for a flat plate, a sphere, and a trans- 
SUBSCRIPT 2 REFERS TO A DIATOMIC GAS verse circular cylinder in a free-molecule flow 
O1 
08 = . ; A 
l 7 ; 1 9y —5 Nu G+hk : 
0.6} a St = = (3.1) 
ayt 1 adt(y + 1) Re 2 
4 A 
fe) y+ 1 l I ; 
03 fmu2 + — = = (3.2) 
¥Y s* G + I 
0.2 This has been done for a flat plate, a sphere, and a 
G transverse flow past a circular cylinder. The algebraic 
and numerical results are represented in tabular form 
hot and are plotted as a function of the molecular speed 
8 . ° . . 
ratio and the free-stream Mach Number for a mon- 
0.06 atomic gas (subscript 1) and a diatomic gas (subscript 
2) 
0.04 ae 3 
oes rable 1 represents the results concerned with the flow 
around the objects under investigation, Table 2 repre- 
002 sents the results referred to their front side only, and 
Table 3 represents those of their rear side with respect 
to the mass flow direction. Fig. 1 represents the plot of 
0.0! the average surface integral G; Fig. 2, that of Ff’; Fig. 3, 





s 02 0304 060810 15 2 3 4 6 8 10 s . . . . 
r ——$—— — —_—$____—— the Stanton Number; Fig. 4, the recovery factor; and 
M, 0.2 03 04 O6 O08 10 15 2 3 4 6 8 10 M os . ‘ . ; ia. 
u ; ee ma ne Fig. 5, the ratio of the adiabatic wall (equilibrium) tem- 
M , a on 
ro Ge Gee WF +5 8 8H perature, 7°,,,, to the free stream temperature, 7). The 
Fic. 1. Average surface integral G as a function of the mo- abscissa of all the diagrams is the molecular speed ratio, 
lecular speed ratio (and the Mach Number) for a flat plate, a ‘ 
sphere, and a transverse circular cylinder. 
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s, with a superimposed scale for the Mach Number, M. 
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Fic. 4. Recovery factor for a flat plate, a sphere, and a trans- 
verse circular cylinder in a free-molecule flow. 
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Fic. 5. Ratio of the adiabatic wall temperature and the free- 
stream temperature for a flat plate: a sphere, and a transverse 
circular cylinder in a free-molecule flow [7 == (Taw/To) — 1). 


The only odd function in the expressions for G and F, 
Eqs. (2.13) and (2.14), is the first term of the F func- 


tion—namely, 


1 { | 
I : dS = 
SJs 2s 2S 


where @ is the local angle of incidence. 
average integral G is therefore the same for the whole 
object as for its front or rear side only; consequently, 
the values of this parameter have not been repeated in 
Tables 2 and 3, and the four curves plotted in Fig. 1 
For the case 


sin 6 dS 
S 


The value of the 


apply to all the cases under consideration. 
of the flat plate parallei to the flow direction, F = 0, the 
heat-transfer characteristics are the same for both sides 
as for each side separately; consequently, this case has 
not been repeated in Tables 2 and 3. Furthermore, 
since the values of the odd term, as shown above, are in- 


dependent of the molecular speed ratio, the value of 
for one side of the object differs only by a constant from 
that for the whole object, for all values of s. 


By virtue 
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of the same argument and from Eq. (2.18), it follows 
that the Stanton Number has the same property. 

The effect of shape on the Stanton Number decreases 
with the decreasing values of the molecular speed ratio 
(or Mach Number), and, for small values of s, the Stan- 
ton Number approaches the value of [(y + 1)/y] + 
(a/4+/ms). Although, for large values of s, the shape 
effect is evident, for high enough values of this parame. 
ter (say s > 2) the variation of the Stanton Number for 
an object of any shape exposed to the free stream cannot 
exceed the value of [(y + 1)/4y] @ [(2/5) @ for a mona- 
tomic and (3/7)a@ for a diatomic gas], and the maximum 
possible local value of the Stanton Number at s = © js 
[(y + 1)/2y]a. 

In contrast to the Stanton Number characteristics, 
the effect of shape on the recovery factor is most pro- 
nounced for small values of s.. For large values of s the 
recovery factor for any surface shape, with the excep- 
tion of the rear side of a flat plate, approaches rapidly 
the value of 2y/(y + 1) (5/4 for a monatomic and 7/6 
for a diatomic gas), so that for s 2 5 it is practically 
equal to this asymptote. The asymptotic value of r for 
the rear side of a flat plate normal to the flow direction 
is zero, while the corresponding asymptotes for the rear 
sides of plates at an angle of incidence lie between zero 
and 2y/(y + 1), as indicated on the right-hand side 
margin of Fig. 4. 

The physical interpretation of this behavior of the re- 
covery factor is not immediately clear. However, since 
the recovery factor serves only for the determination of 
the adiabatic wall temperature, 7°,,,, and is used only in 
conjunction with the Mach Number or the molecular 
speed ratio, the situation becomes clearer if one con- 
siders the ultimate result of the recovery factor applica- 
tion by observing the behavior of parameter 7, defined 
as 

T ase Y , 


t=—— 1 = ie (3.3) 
To Y 

which is plotted on Fig. 5. The ordinate of this diagram 
is expressed as [(y + 1)/(y — 1)] 7 so that its value is 
simply equal to the product of corresponding ordinate 
of Fig. 4 and s? {since [(y + 1)/(y—)]r = [(y #1) + 
y]rs?}. As evidenced by Fig. 5 and as should be ex- 
pected, the effect of shape on the ratio of adiabatic wall 
temperature to free-stream temperature is small for low, 
as well as for high, values of the molecular speed ratio 
and is most pronounced in the vicinity of s = 1. 


(3.2) Applications 

Besides the shapes analyzed so far and indicated on 
the diagrams, the heat-transfer characteristics of a flat 
plate inclined at any angle of incidence, 6, to the flow 
direction can be easily found by extrapolating the results 
corresponding to the flat plate normal to flow, since, by 
Egs. (2.13) and (2.14), respectively, 








G(s) 


and, $ 


where 
paratl 
incide 
plate 

Fre 
(3.3), 


On 
prese 
trans 
incid 
the f 


since 


F 
sult: 
cide 
flat 

thos 
cont 
trar 
are 

any 


Syn 





lows 





reases 
ratio 
Stan- 
y] + 
shape 
rame- 
er for 
annot 
10na- 
mum 


© 15 


stics, 
pro- 
5 the 
Cep- 
idly 
7/6 
‘ally 
for 
tion 
rear 
eTO 


side 


nce 
of 








CONVECTIVE HEAT TRANSFER 
id ne” 
Sil we oS) = - -=- =- = sin 0G,(n) 
aa? * SH" Q/ 4s s 2 xn — 
(3.4) 
and, similarly, 
F,(s) = fo(s) = sin 6 F,(n) (3:5) 


where » = S sin 0, subscript @ denotes the value of the 
parameter corresponding to the flat plate at an angle of 
incidence @, and subscript » denotes that of the flat 


plate normal to flow direction. 


From the above by virtue of Eqs. (3.1), (3.2), and 
(3.3), it follows that 
Ste(s) = sin 6 St,(n) (3.6) 
9 
r(s) = , cos” @ + sin? 6 7,(n) (3.7) 
Y 
mee ; 
T)(Ss) = : s? cos? 6 + 7,(n) (3.8) 


On the basis of these expressions and using the results 
presented here in tabular or graphical form, the heat- 
transfer characteristics for a flat plate at any angle of 
As an example for 


the front side of a flat plate at @ = 30° and s = 2, 


incidence can easily be computed. 


St = 0.5 X 0.5126 = 0.2563 
Qa y + l 
a | a 
eT pw id + 0.25 X 2.8988 = 2.2247 
y 
since, from Table 2, forn = 2 X 0.5 = 1, 
I ef ; ete 
St,(1) = 0.5126 
aytl 
I= r, (1) = 2.8988 
y 


Furthermore, it should be noted that the flat plate re- 
sults apply to any surface with a constant angle of in- 
cidence—that is, the heat-transfer characteristics of a 
flat plate at an angle of incidence @ are the same as 
those corresponding to the axial flow past a wedge or a 
cone with @ as the semiapex angle—and the heat- 
transfer data of a flat plate tangent to the flow direction 
are the same as those for axial flow past a cylinder of 


any shape exposed to the free-molecular stream. 


TABLE 4 
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Fic. 6. Illustration for numerical example. 


Finally, the results of the analysis can also be extra- 
polated to other surfaces by noting that, by virtue of 
Eqs. (2.13), (2.14), and (3.1), 


St = (1/S)>>St,S; 


¥+1 l SF,S, 
r=2+ . ao 
7 s GS; +> DF;S; 


where subscript 7 denotes the fundamental components 


(3.9) 


(3.10) 


of a given surface. 

As an illustration, the preceding example is con- 
tinued in order to find the heat-transfer characteristics 
of a simple body shown on Fig. 6 at s = 2. Basic values 
for computations are presented in Table 4, and the re- 


sults are as follows: 


l 7 0.1987 5 
= ———- «= 6258 
Pars 1.75 
r+] | 0.22975 
: —s rina te SEO 
. 2? (0.1676 + 0.22975 


SUMMARY 


The analysis of the convective heat transfer in a free- 
molecule flow with Maxwellian velocity distribution has 
been generalized by reducing all the relationships to a 
form common to all gases irrespective of their molecular 
structure. The heat-transfer characteristics of a num- 
ber of simple surfaces consisting of a flat plate in the 
direction of flow, a flat plate normal to flow, a trans- 
verse cylinder, and a sphere have been computed and 
method has been outlined for 


correlated. A simple 





l 2 3 4 5 6 7 

S/r)D? G F St S-G S:F S:St 

Symbol Source: Fig. 6 Table 1 Tables Table 3, pre- 1-2 1-3 1-4 

2,1,3 vious example 

Si Front cone 0.5 0.0519 0.4607 0.2563 0.02595 0.23035 0.1282 
Se Cylinder l 0.1410 0 0.0705 0.1410 0 0.0705 
S3 Back plate 0.25 0.0026 —0.0024 0.0001 0.00065 — 0.0006 0.0000 
0.22975 0.1987 


tw 
“J 
or 


0.1676 
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extrapolation of the results to a flat plate at any angle 
of incidence, a wedge, or a cone and for the evaluation of 
heat transfer parameters for a compound surface. 

The effect of shape on convective heat-transfer char- 
acteristics has been found to be small. At low Mach 
Numbers, the Stanton all shapes ap- 
proaches the value of 


Number of 


[vy + 1)/y]-(a/4-/ xs) 


The difference between various shapes becomes more 
distinct at high Mach Numbers. For J/ > 2, the mag- 
nitude of the Stanton Number cannot exceed the value 
of [(y + 1)/2y]a. The effect of shape on the recovery 
factor is more distinct at low Mach Numbers and be- 
comes hardly noticeable for Mach Numbers higher 
than 5. For the temperature ratio, 7\,,,/7», this effect 
is small for both low and high Mach Numbers and is 
most noticeable in the vicinity of Mach Number 
unity. 
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RIEF REPORTS of investigations in the aeronautical sciences and discussions of papers pub- 


lished in the JOURNAL are presented in this special department. 
The 
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for the opinions expressed by the correspondents. 


A Rational Discussion of the Hydraulic Analogy 


E. V. Laitone 
Associate Professor, University of California at Berkeley 


September 5, 1952 


T= NOTE ISAN ANSWER toreference 1,acomment on my paper- 

on the application of the hydraulic analogy to transonic flow 
studies. Reference 1 doubts my statement that the analogy is 
restricted to '/,-in. water depths since the group velocity, which 


is given by*® 


— Ndc/dd) (1) 


a= tC 


approaches the wave velocity c as \ ~ ©. However, who can 


build an infinite model, much less move it ? 


For the analogy to hold, the group velocity must be equal to 
the wave velocity, as shown in reference 2. Now for ordinary 
> a for all h > '/, in., 
the 


finite wave lengths A, we have « as may be 


shown by neglecting the surface tension. Then so-called 


wave velocity (which is actually only a phase, net a true propaga 


tion, velocity ) is given approximately by 
P I/» ‘ 
c = [(g\/27) tanh (2rh/d)] (2 


Then, combining Eqs. (1) and (2) we have the ratio of the group 


velocity (a) to the wave velocity (c) given by 


a r (* l ( toh/d 
nk Fa - = (3) 
c c \dy 2 sinh 4rh/X 


so that 
for h/X >> 1 and 


for h/A< 1. 


10h have a varying group velocity 
For ex 


This proves that all A 
ranging from 50 to 80 per cent of the wave velocity. 
ample, a water depth of 1 in. requires that all the wave lengths 
less than 1 ft. be practically eliminated by using a sufficiently large 
model. However, making a finite wedge, no matter how sharp 
the nose angle, 1 ft. wide does not mean that only wave lengths 
long have a significant influence on the flow pattern, 
Consequently, the model would 


over 1 ft 
especially in the transonic case. 
have to be enormous in size before the analogy became valid 


Comparing the effect of water depths over ‘/, in. for the case 
of an attached oblique shock wave at the nose of a wedge does 
not even enter into a discussion of the effect of group velocity. 
This finite discontinuity is a hydraulic jump having a momentum 
discontinuity that has no direct relation to the group velocity de 


fined by Eqs. (1) and (2). The group velocity is a potential-flow 


Publication will be completed as 
Editorial Committee does not hold itself responsible 
Contributions should not exceed 800 words in length 


relationship that gives the analogy to a potential gas flow with 
y = 2, while the hydraulic jump is an irreversible (nonpotential ) 
finite discontinuity that does not have a direct analogy with a 
shock wave in a perfect gas. However, weak shock waves ap- 
proach the potential-flow solution; consequently, the analogy is 
valid for weak discontinuities as occur in transonic flow or for 
very slender wedges. This again means that the water depth 
must be of the order of '/, in. 

It is interesting to note that reference 4, the same report as re- 
ferred to in reference 1, shows that even the experimental data 
corresponding to finite oblique shocks are in much better agree- 


ment with a perfect gas having y = 2 when the water depth is 
about '/; in. (Fig. 8, reference 4). Asa matter of fact, the author 


states on page 19 of reference 4 that the oblique shock angles 
measured in water compared most favorably with that for a per 


fect gas having y = 2 when the water depth was 0.225 in 
& g 


REFERENCES 


| Randels, W. C., Comments on ‘‘A Study of Transonic Gas Dynamics by 


the Hydraulic Analogy,’’ Readers’ Forum 
19, No. 8, p. 572, August, 1952 
Transoni 


Journal of the Aeronautical Sci 


ences, Vol 

? Laitone, E. V., A 
{nalogy, Journal of the Aeronautical Sciences, Vol. 19, No. 4, pp 
April, 1952 

* Lamb, H 
England, 1932 
4 Bruman, J. R., 
ey, North American Aviation, Inc., Report No 


A ydraulic 


265-272 


Study of Gas Dynamics by the 


Hydrodynamics, 6th Ed., p. 381; University Press, Cambridge, 


A pplication of the Water Channel-Compressible Gas Anal 
NA-47-87, March 3, 1947 


A Heat Cycle Low Reynolds Number Transonic 
Tunnel 


Raymond L. Chuan 
Graduate Student, Guggenheim Aeronautical Laboratory, 
California Institute of Technology, Pasadena, Calif 


August 11, 1952 


Pp eeretey WIND TUNNEL with low Reynolds Number flow 
designed to aid the study of viscous effects has been built 
at GALCIT. 
operates at Mach Numbers up to 1.8 and at Reynolds Numbers 


Using water vapor as the working fluid, the tunnel 
of about 5,000 per cm. The continuous operating closed return 
circuit tunnel has no mechanical compressor 
vapor in the diffuser is condensed; the liquid flows under the in 
fluence of gravity to the stagnation chamber where it is evapo 
The pressure differential between diffuser 


The superheated 


rated and superheated. 
and stagnation chamber ordinarily supplied by a mechanical com 
pressor is maintained in the new tunnel by a simple difference in 
liquid level between these two chambers. The low stagnation 
chamber pressure (20 to 35 mm. Hg) and density (10~® Gm. per 
the Reynolds Number 


cu.cm.) necessary to produce desired 
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would have necessitated operation at extremely low efficiency 
had a mechanical compressor been used. 

For a tunnel with fixed Mach Number, test-section dimensions, 
and stagnation temperature, the mass flow decreases with the 
Reynolds Number, 
Hence, if mechanical compressors are used, the physical size of 


while the volume flow remains constant. 
the compressors remains about constant, and the power con- 
sumption of the compressors decreases more slowly than the mass 
flow. 
change system, as in the new tunnel, the power consumption 


On the other hand, if the compression is done by a phase- 
decreases with the mass flow. One advantage of a phase-change 
compression system for a low Reynolds Number transonic wind 
tunnel is then obvious: low power consumption. 

This principle of using phase changes in the compression proc- 
ess was suggested to the author by Prof. H. W. Liepmann in the 
summer of 1950. 
ingly been built, and the initial success indicates that the feasi- 


The present water-vapor tunnel has accord- 


bility of wider application of the principle warrants further study. 
For instance, the application of the phase-change compression 
principle to hypersonic tunnels appears especially attractive since 
the extremely high compression ratio necessary for hypersonic 
Mach Numbers can be accomplished at nearly 100 per cent effi- 
ciency with leakproof, unlubricated liquid pumps. 

Fig. 1 shows the complete tunnel. The stagnation pressure 
and temperature are regulated by two autotransformers con- 
trolling the power input to the electrical immersion heaters and 
superheaters. The nozzle has flexible walls controlled by two 
pairs of jacks and provisions for model mounts and traversing 
The condenser is driven by a refrigerator that re- 
jects heat into cooling water, although power consumption could 
In 


addition to the conventional pressure switches on the refriger- 


mechanisms. 
be reduced by rejecting part of the heat into the boiler. 


ator, a by-pass evaporator is provided for regulating the heat re- 
moval in the condenser. The difference in water levels between 
the condenser and the boiler makes up the pressure difference be- 
tween the diffuser and the stagnation chamber. 

Instrumentation of the fluid 


manometers to measure pressure referred to the stagnation pres- 


consists following: silicone 


sure; absolute pressure gage to measure stagnation pressure; 
dial-type thermometers to measure temperatures of the boiler, 
superheater, and condenser; and gage glasses for observing liquid 
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The low Reynolds Number tunnel. 
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Fic. 2. Flow at M = 1.7 over a 20° wedge. 
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Fic. 3. Mach Number profiles in test section 
levels. A two-mirror schlieren system, arranged with symmetry 


shock 


*Gm./cm.‘, which are about 


to minimze distortions, is used to detect waves with 
density gradients of the order of 10 
five orders of magnitude below those in corresponding shock 


waves in more conventional supersonic tunnels. 


Typical operating data are as follows: 


po = stagnation pressure = 30 mm. Hg 

Ty = stagnation temperature = 180°C. 

po = stagnation density = 2 X 10-5 Gm. per cu.cm 
p = static pressure in test section = 6 mm. Hg 

M = Mach Number = 1.7 

p = static density = 5 X 10-6 Gm. per cu.cm. 

A = mean free path ~ 1073 cm. 

R,. = Reynolds Number = 4,500 per cm. 


W = total power consumption = 7 kw. 
A = test-section cross-sectional area = 3.9 sq.cm. 
Fig. 2 shows a schlieren photograph of the flow over a 20° 


ee 
across the test section. 


wedge at M = Fig. 3 shows two Mach Number profiles 


At present, two problems are being investigated using this 
tunnel: One is the curvature of an attached shock wave in the 
close neighborhood of the nose of a wedge, with the iutention of 
deducing some knowledge of the boundary-layer structure there; 
the other is a study of heat transfer from a heated flat ribbon, 
which would be difficult to carry out in a tunnel at higher density 
because of the greater possibility of vibration difficulty at high 
dynamic pressure. 


The initial success of this low-density tunnel has not only pro- 
vided an instrument for the study of flow in the lower part of the 
continuum gas dynamics régime but has pointed the way to wider 
applications of the vapor cycle to other wind tunnels. Further 


study of the latter aspect is proceeding. 
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On Second-Order Supersonic Flow Past a Slowly 
Oscillating Airfoil 


Milton D. Van Dyke 
Ames Aeronautical Laboratory, N.A.C.A., Moffett Field, Calif. 
October 15, 1952 

HE WRITER HAS RECENTLY CALCULATED the second-order 
Tae of thickness upon a slowly oscillating airfoil in super 
sonic flow. Because the results disagree with those of other 
treatments of the problem, » ? it seems worth while to summarize 
them here 
Consider a sharp-nosed airfoil whose upper surface is described 


at zero angle of attack by y = V(x). Let the airfoil pivot about 
a point lying a distance ) downstream from the leading edge and 
oscillate with small angle of attack a(t) (for example, a = ap sin 


wf for sinusoidal oscillations). A solution in terms of elementary 
functions is obtained by restricting consideration to slow oscilla 
tions so that, for example, in the case of sinusoidal motion only 
linear terms in frequency are retained. Then, to second order in 
thickness (but only first order in angle of attack, which is suffi 
cient for investigating stability), the pressure coefficient on the 


upper surface is found to be 


P. 2a (2 — M? 
=-(Y’—a)+ ; x+b)+ 
8 Bl B? 


C, = 
WwN—2 : : 2a j N-1 
(yY’? - 2} a) + J 4 2M? } T 
B2 4 4 

(2 — M*)(AMP2N — 1) M?@N —2_ | 
x 0 r (1) 

B4 B? f 
Here, U and M are free-stream speed and Mach Number, 6 = 
(M? — 1)'/2, and N = (1/2)(y + 1)M?/6?, where y is the adia 


batic exponent. The prime and dot indicate differentiation with 
respect to x and time ¢, respectively 

The first two terms in Eq. (1) are the usual result of linearized 
theory. The third the 
term, which would result from using the instantaneous slope in 
The fourth term is the new result of this 


term is second-order ‘‘quasi-steady”’ 
Busemann’s formula 
analysis 

Consider now the stability of pitching oscillations for a sym 
metrical airfoil with sharp trailing edge. The pressure on the 
lower surface is obtained by reversing the sign of @ in Eq. (1 
Thus, the pitching moment coefficient (about the pivot) is found 
to contain, in addition to a restoring term, a damping term given 


by 

4ca}2— MM? 3b bh b M%7N — 1) 

— _ } 1—2-})- . x 

I 383 pe 281 ( B' 

i awn -2)1 {* 

- Y dx +4 (x — b)Y dx (2 
 J0 p* c* Jo 

where ¢ is the airfoil chord. The moment is stabilizing or 


destabilizing according as this term is negative or positive. It is 
well known that linearized theory, represented by the first two 
terms in brackets, indicates the possibility of instability for cer- 
tain combinations of Mach Numbers less than (5/2)'/? and pivot 
points ahead of two-thirds of the chord. The effect of thickness 
is to shift this range of possible instability toward higher Mach 
Numbers and more forward pivot points, as shown in Fig. 1 for 
a6 per cent thick double-wedge airfoil or 4 '/2 per cent thick 
biconvex airfoil. (The dotted portion of the curve corresponds 
to subsonic flow behind the bow wave for either airfoil, so that 
the theory is inapplicable. ) 

The second-order effects of thickness are small, as would be 
expected, and ordinarily tend to destabilize. These results differ 
from those of reference 1 and are in sharp disagreement with the 








1eT 


results of reference 2, which predict strong stabilizing effects 
of thickness. The details of the theory of reference 1 are not 
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available for discussion. The analysis of reference 2 appears to 


contain a basic error. The author avoided the difficulties associ 
ated with kinked streamlines by the strategem of considering a 
general ‘“‘smoothed”’ problem and then passing to the limit of 
sharp corners in the solution. This artifice is justified because 
to second-order a shock wave is the limit of a rapid continuous 
isentropic compression. However, it appears. that the author 
failed to solve the actual smoothed problem before passing to the 
limit. 

Just at the 


nose of an oscillating airfoil (so long as the shock wave remains 


A simple test can be applied to these solutions. 


attached ), the pressure depends only upon the instantaneous local 
velocity (see reference 3). Therefore, at the nose the pressure 
must equal that on the same airfoil in steady flow with the in- 
stantaneous angle of attack increased by the apparent downwash 
angle, ab/U. The [Eq. (1)] this test, 
whereas the result of reference 2 fails to do so. 


present result passes 


The writer expects to publish soon the details of his analysis. 
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A Method for Location of Detached Shock 
Waves Ahead of Plane Bodies 


Robert Wesley Truitt 

Professor of Aeronautical Engineering, Virginia Polytechnic 
Institute, Blacksburg, Va. 

September 8, 1952 


DISCUSSION 


ECENT EXPERIMENTS ON WEDGE PROFILES, reported in refer- 

ence 1, reveal that in transonic flow the shock detachment 
distance is affected considerably by the opening angle of the pro- 
file at the leading edge. It is the purpose of this paper to derive 
an expression for the shock detachment distance which is an ex- 
plicit function of the opening angle of the wedge 
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Fic. 1. 


The pertinent limiting conditions for transonic flow are par- 
ticularly well defined for the present study and occur at JJ) = 1 
and at minimum shock attachment Mach Number, J/g, for a 
given profile. Considering the fact that sonic velocity occurs at 
the shoulder of a wedge (see reference 2) and using the limiting 
conditions for shock detachment at 1/) = 1 and J/g, the shock 
detachment distance, S, for any transonic Mach Number may be 
found (see Fig. 1). 

Consider the angle 6’ in Fig. 1, where the first limiting condition 
The other limiting 
An ex- 


is such that 0’ > 0 as Ip 1 (i.e., S> @), 
condition is such that 0’ ~ @as My —~ Mg (i.e., S—> 0). 
pression that satisfies these limiting conditions may be written as 
(see Fig. 1). 

tan @’ = t/(S+L)Ax0’ 


or 
0’ = (t/L)/[(S/L) + 1] (1) 

where 0’ may now be identified as the angle corresponding to the 
minimum free-stream Mach Number, J, (1 < AJ) < Mg), which 
will still produce an attached shock wave at the nose of a wedge 
having a nose angle of 26’. In other words, exact shock theory 
gives for every transonic Mach Number, Wo, the maximum nose 
angle 6’ which will still allow shock attachment. 

From Fig. 1, {/L & 6 (the semiangle of the wedge), and Eq. (1) 
may be written in terms of the shock detachment parameter, 
S/L, as 


S/L = (6/0’) — 1 (2) 
Again, when My = Mg, 6’ = @and S = 0; when M = 1, 


6’ = OQand S = ». From the above consideration it appears that 
Eq. (2) should give reasonably accurate values of S/Z in the en- 


tire transonic range. 


“i — 


2.0+ 


---0-- THEORY- EQ.(2) 
————__ EXPERIMENT 
(Ref. 1) 





S/L 
12+ 
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COMPARISON WITH EXPERIMENT 


Eq. (2) has been used to calculate S/L values in the transonjc 
range for 41/2°, 71/.° 
The calculated values of S/L have been compared with the ex. 


perimental results of reference 1 and are presented in Fig. 2 


and 10° two-dimensional semiangle wedges 


The agreement of Eq. (2) with these experiments is evident 


CONCLUSIONS 

Other simple expressions have been derived for shock detach. 
ment distance which predict no direct effect due to the nose shape 
or varying nose angle.* This prediction is not true, as is borne 
out by experimental evidence, for thin wedges where the detached 
shock location is considerably affected by the nose angle Eq 
(2) predicts a considerable effect due to nose angle, as is evident 
from Fig. 2. The present analysis has been made only for the 
two-dimensional case of thin wedges with well-defined shoulders 
and had not been attempted for bodies of revolution at the time 


of writing. 


REFERENCES 


E., Jr., An Experimental Investigation of Transonic Flow Past 
Two-Dimensional Wedge VM ach-Zehnde 
Interferometer, N.A.C.A. T.N. No. 2560, November, 1951 

2 Cole, J. D., Drag of a Finite Wedge at High Subsonic Speed 
and Phys., Vol. 30, No. 2, July, 1951 

3 Moeckel, W. E , Approximate Method For Predicting Form and Location 
of Detached Shock Waves Ahead of Plane or Axial 
N.A.C.A. T.N. No. 1921, July, 1949 


' Bryson, A 
and Circular-Arc Sections Using a 


, Jour. Math 


Symmetric Bodie 


—— 


A Note on Tabular Methods for Flexural 
Vibrations 


William T. Thomson 
Department of Engineering, University of California at Los Angeles 
September 24, 19592 


= METHODS OFFER a practical approach to problems of 
beam vibrations where the stiffness and mass distribution are 
nonuniform. These methods are all based on the fact that, in 
forced harmonic oscillations, the deflection, slope, moment, and 
shear are linearly related along the beam. It appears that such 
problems can be formulated in many different ways, and the pres- 
entation here offers an advantage of algebraic simplicity which 
retains the physical concepts of the beam. 


UNCOUPLED FLEXURAL VIBRATION 


The beam section is chosen as in Fig. 1. The influence co- 
efficients are defined differently from those of previous papers as 
follows: The coefficients as, ay, ds, and dy are the slope or de- 
flection at station 7 +4 


toa unit shear or moment at 7. 


1 measured from a tangent at station 7 due 














Fic. 1. Beam section for uncoupled flexural vibrations. 
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With the influence coefficients so defined, the computational 
ons for the uncoupied flexural vibrations can be obtained 


equati 

directly from Fig. l as 
Sign = Si + miw?y; (1) 
Mins = M; + Sisal; (2) 
Gin = Gi — Sizqiasi Miami (3) 
Virol = 9; ¢ili + Siuidsi + Miday: (4) 


where the first two equations are deduced from equilibrium con- 
siderations and the last two are from the geometric configuration. 

For a chosen value of w, these equations enable one to calculate 
5, M, ¢, and y at any point in terms of the corresponding quanti- 
ties at some starting point. For instance, for a free-ended beam, 
the computations are started with the following boundary values 
at the free end: 


S=Mm=0, ga=¢, wm = 1.0 


Progressing to the other end, the quantities obtained are in the 


form 
Se = a, + bigy (5) 
Mn = a2 + bog, (6) 
gn = 03 + b3¢1 (7) 
Vn = 4+ bay (8) 


where the a’s and b’s are numbers. Specifying the boundary con- 
ditions at ”, ¢: is determined, after which all other quantities at 
n are established. Procedures for the various boundary con 
ditions are well known and are not repeated here. 

Eqs. (5) to (8) indicate that the a’s and b’s independently satisfy 
the computational equations, Eqs. (1) to (4), thereby suggesting 


atabular form in two parts 


ROTATING BEAMS 


The forces acting on a section of a rotating beam are shown 
in Fig. 2, where the vibration is assumed to take place in a direc 
The computational 


established from 


tion perpendicular to the plane of rotation. 
equations, including the effect of rotation are 
this diagram to be 


Sinn = Si + miw*y, — mi Q?rigi (9) 
Migs = Mi + (Si + myw?y;)l; mQ?r,(y¥;5 — Via) 

= Mil + m22ridyi) + Sih, + mQ?rids;) (10) 
Gin = Gi Siqias Miami (11) 
Hn = Vi—- © lL, + Sjiuids: 4 Mid Mu (12) 


It is evident here that only the shear and moment equations are 
altered in form by rotation and that the computational procedure 


and tabular form for the uncoupled flexural case apply. 


CouPLED FLEXURE-TORSION VIBRATION 


Natural modes of vibration of airplane wings are always coupled 
flexure-torsion modes, which for higher modes may differ con- 
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Fic. 2. Beam section under rotation. 
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Beam section for coupled flexure-torsion vibrations. 
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Fic. 4. Stiffness curve for determination of influence coefficients. 


siderably from those of uncoupled modes. It is necessary here 
to introduce one additional influence coefficient h;, defined as the 
angle of twist at station 7 + 1 with respect to station 7 due toa 
unit torque at station 7. Referring to the beam section shown 


in Fig. 3, the following equations can be written: 


Sinn = Si + miw*(yi + €:6;) (13) 
Mia = Mz + Sig ls (14) 
Tia. = Ty + Jw; + m Cw"; (15) 
Yi+l = Yi Si41 Qsi Miami (16) 
Vint = VM gili + Sins dsi + Midas: (17) 
Gin. = 0: + Tih; (18) 


For free-ended beams we have the following boundary condi- 
tions to start with: 


g_a=¢, N= 10, W 


= 0 


Here, again, the quantities of interest at any station are linearly 
related to yg; and 6; and can be expressed in the form 


at be, + 6; 


This suggests a tabular form where each quantity 5, /, 7, etc., 
would be divided into three parts—namely, a, 6, and ¢. 
INFLUENCE COEFFICIENTS 


The influence coefficients as, ay, ds, and dy are determined 
from the actual stiffness curve of the beam by the area-moment 
procedure. Letting the ordinates at i and i + 1 be aEip and 
(a + b)EIy, where El is a reference stiffness, the curve between 
the two stations is replaced by a straight line, which results in the 
following equations: 

a I “ dé 


0 El El, Jo (a+b) — (b/E 


1; (° + 
In 
hEI a 


ami = 


b 
) (19 
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“~~ ee 1,2 a+b (2em + fe) fee + 2ca fer — Can.fan — fre = 0 (4) 
as; = = = = : 6 —-— ah} —— (20) ; ; ; : ; 
0 El BEI, for the dimensionless perturbation potential /. Previous special 
Li cases of this equation have been given by Oswatitsch and Wies. 
& dé 1;2 a+b a q : _ Ore d Wieg 
dui = - = - (a + b) In —}—>p (21) whardt? (cy = cy = Ca.r. = 0), von Karman* (cy, = Ca. rn. = 0), ang 
0 El BEL a Gat texa 0). 
“i (1; — te dé The similarity parameter c,.p, was introduced first by Berndt,' 
dsi = - oe —— = but he implied that the potential equation may be linearized only 
0 ‘ in the limit A.R. = 0, whereas ca.p. can be large for only moder. 
1,3 l a+b : itil cidin =e ; = ae = pee. 
2 oY + ob — cla + 0a (22) ately small values of A.R. (For 6 = 0.05 and A.R. = 1, we haye 
b3EIo L2 a Car. = 4.1.) If car. is large and both cy and c, are 0(1), we 


It is of interest to compare the influence coefficient as defined 
here to those of Myklestad 
measured from a tangent at 7 + 1 due toa shear or moment at 7. 
Using the bar notation for the latter, the relationship between the 


i.e., slope and deflection at station 7 


two are 


ami = &yi 
ag, = &35; = dai 
dy; = ds; — leu; 


dsi = dsi — lidmi 


The above relationships indicate that the influence coefficients 
as used here can be checked experimentally in the usual manner. 
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On Linearized Theory 


John W. Miles 
Department of Engineering, University of California at Los Angeles 
September 8, 1952 


HE NECESSARY CONDITIONS for the linearization of the equa- 
tions of irrotational compressible flow past a thin two- 
dimensional body have been set forth by Lin, Reissner, and 
Tsien.! 
sional wings of thickness ratio 6 and aspect ratio A.R. shows that 


An extension of their analysis to flow past three-dimen- 


necessary conditions for linearization in the neighborhood of the 
wing are that all of the conditions 


6, Mb, kb, RM5< 1 (1) 
be satisfied and that any one or more of the conditions 
M-1|>6" (2a) 
k>s’s (2b) 
1/A.R. > 6/3 (2c) 


be satisfied. M is the free-stream Mach Number, and & is the 
reduced frequency (or, more generally, a dimensionless measure 
of time ratio of change). In the case of a slender body of revolu- 
tion,? we may set A.R. = 6, sothat Eq. (2c) is included in Eq. (1). 

Conditions equivalent to Eq. (2) may be obtained by introduc- 


ing the transonic similarity parameters 


cw = [(y + 1)8]7~*7|M —- 1 (3a) 
ce = [((y + 1)6]~ (3b) 
CAR. = [(y oo 1)5] + (A.R.)-? (3c) 


and requiring one or more to be large (relative to unity). On the 
other hand, if ail three parameters are 0(1) we obtain the dimen- 


sionless nonlinear transonic equation [cf. Eqs. (21) and (25), 


reference | | 


have Jones’ slender wing theory.’ It should be remarked that 
Ca.R, >> 1 is a necessary condition for slender wing theory in the 
transonic régime, whereas JJ = 1 (socy = 0) is neither necessary 
Similarly, while (17? — 1 )'/2A_R. is a valid tran- 
sonic similarity parameter |[=(cy/ca.r.)'/*], it is neither necessary 


nor sufficient. 


nor sufficient that it be small for the linearization of Eq. (4). 

It is not possible to assert that the foregoing conditions are 
sufficient without (a) imposing restrictions on the curvature of 
the wing surface and (b) investigating the boundary conditions 
at infinity. In general, the linearized approximation will not be 
uniformly valid in the neighborhood of the wing edges, but this 
difficulty can be circumvented by the method of Poincaré and 
Light hill.’ 
disturbance of finite energy, and in transonic flow past an airfoil 


At infinity the solution must behave as an outgoing 


at zero incidence, this can be ensured only by requiring the asymp- 
totic solution to satisfy Eq. (4), even though Eq. (2c) may per- 
mit linearization in the neighborhood of the wing 
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Interaction Between Stalled and Unstalled Flow 


Peter F. Jordan 
Now Royal Aircraft Establishment, Farnborough, Hants., England 
October 1, 1952 


— IN A PAPER recently referred to by Bird? discussing ex- 
perimental lift distributions over a A-wing in partly stalled 
flow, writes: “The spanwise variation Of Cangz, is not explait- 
able... . since Camar. reached a maximum at 39.6 per cent semi- 
span rather than at the center of the wing .. . It can be con- 
cluded that there exist as yet undefined effects of the boundary- 
layer flow on the wing stalling characteristics.” 

The writer? investigated the mutual influence of stalled and 
unstalled flow (AVA-Géttingen, 1937 and early 1938). He found 
that excessive sectional cj-values were due to periodical occuf- 
rences. Let y be the spanwise coordinate and let yo be the border 
between stalled and unstalled flow; let y > yo be the stalled region. 
Then the flow will suddenly attach itself at a point y, > yo, thus 
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forming an island of unstalled flow; this island will immediately 
move toward yo, making room for a new island to be formed at ,, 
and soon. The start of this process is explained partly by a 
minimum of induced downwash that exists near y, and partly by 
a spanwise boundary-layer flow near L.E.; the process as a whole 
produces, between yo and y,, an effective camber of the wing 
profile, thus leading to excessive c;-values. 

The above observations were made on wings with constant 
chord; of course they are not sufficient to explain fully the more 
complicated occurrences on A-wings, but they seem to explain 
why Cnmar. Will occur neither at mid-span nor at the wing tips 
Indeed, the c;-value belonging to Cp,,,,, as measured by Wick! 
was far in excess of the two-dimensional ¢;,,,,-value for his wing 
profile in two-dimensional flow (1.7 against 0.9). Such excessive 
lift coefficients require interaction between unstalled flow and 
stalled flow. 

It should perhaps be mentioned that a boundary-layer fence 


for the first time to the writer’s knowledge—to de- 


was used 
termine the part played by the boundary-layer flow; it 
shown to effectively fence in the stalled region. The Reynoids 
10°, a rather low value; 


was 


Number of the test was Re = 0.2 X 
however, the flow observations were repeated and confirmed on an 


aircraft in flight by Liebe.* 
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Physical Explanation of Singular Integrals at 
Supersonic Speeds 


Lee S. Wasserman 

September 19, 1959 

Chief, Aerodynamics Research Branch, Flight Research Laboratory, 
Wright Air Development Center, Wright-Patterson Air Force 
Base, Ohio 


HIS NOTE PRESENTS a physical explanation of a type of singu 
lar integral appearing in linearized supersonic theory. We 
consider a pressure source at the origin expressed by Eq. (1) 

(1) 


(x,y,z) = €/2n V x? — (M? — 1) (y? + 2?) 


The vertical acceleration of an air particle at x, y, z, due to the 


pressure source is given by Newton’s Law. 


dw Op eM? — 1)z 
p —_—<— = — (2) 
dt Oz 2r[x? — (M2? — 1)(y? + 2?)] 
For steady flow dW/dt = U(dw/dx), so that integration yields 
for w 
€ . XZ 
w= - (3) 
2rpU | 2x (y? + 3?) [x? — (M*? — 1) (y? + 2?)] ”*) 


Where x* is the value of x which makes the denominator vanish 
so that the velocity is infinite. 

The trouble here comes from the fact that the derivative taken 
in Eq. (2) does not include the positive pressure jump that occurs 
just ahead of the Mach cone. It is possible to show that the ve- 
locity induced by this jump cancels out the infinity in Eq. (3) 
From acoustic theory, the velocity normal to a plane sound wave 


is given by the expression p/pc. The wave front on the Mach 
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cone appears plane to particles very close, so that w due to the 


jump is given by the expression 
p z cos B 
pc V y? + <2 


ez cos B 
9 (4) 
onpe v/y* + 3s V/x*? — (M2 - 1) (y? + 2?) 


where z cos B/V y? + 2? gives the component of the axially sym- 
metric induced velocity in the vertical direction. Now 
—_ : 
cos B/c = VY M? — 1/1 
/ / 
x*/V y2 + 3%U) = V M* -— 1/U 
so that the velocity induced by the pressure jump in Eq. (4) ex- 
actly cancels the infinite velocity given by Eq. (3). 


—+ 


Frequency Response and the Transfer 
Functions of the Human Pilot Measured in 
Flight 


G. M. Andrew 
Air Force Flight Test Center, Edwards, Calif. 
August 11, 1952 


INTRODUCTION 


“— IMPORTANCE OF determining experimentally the transfer 
function of the human pilot was recognized only in the re- 
cent years. Some attempts of measurements along the line were 
made.' 2 In these tests flight simulators were used. 

The method proposed in this note is based on the actual flight 
tests. The pilot reacts thus under actual flying conditions, and 
such factors as altitude, accelerations, fatigue, speed, etc., are not 
neglected. 

The theoretical background of the method, based on the theory 
of servomechanism, is given in detail in reference 3 

The transfer functions of the human pilot should be determined 
about all three axes. For purpose of illustration, only pitching 
axis is discussed below. 

INSTRUMENTS AND FLYING TECHNIQUE 

A simple angle of attack indicator can be used if gyro equipment 
is not available. 

The time history of angle of attack a must be recorded 
an oscillograph or some other recording device is required. 
time history of elevator deflection 6, must also be recorded. 

The method requires little flying time. Two flight-test records 
are needed for each configuration. 


(1) During the first flight test the pilot puts a step or, still 
The recorded elevator 


Thus, 
The 


better, a single pulse elevator motion 
motion and the aircraft response are shown in Fig. 1. 


f 





6, 


trim 


rim 





FIG. | 
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(2) The panel instrument used by the pilot for observation of 
the airplane response should have a visible mark for the trim 
position of the angle of attack a trim. Immediately after the 
first flight test is finished, the second flight test can be started. 
Different flying techniques can be used, such as starting from 
initial steady state and trying to reach another selected steady 
state marked on the instrument. 

Or a dynamic method can be used, as one described below. 
The pilot puts a considerable elevator step disturbance (Fig. 2). 
At time tz, when the airplane response reaches its maximum value 
ay, the pilot intervenes with elevator motion, trying to get back 
as fast as possible to the marked trim position atrim. The motion 
of the elevator and of the airplane response are, of course, re 
Both flight tests (for one configuration) take only few 

It is important to tabulate gross weight, c.g. position, 


corded. 
minutes. 
altitude, speed, and other pertinent data describing the config- 
uration. 
DATA REDUCING 

Fig. 3 represents the block diagram of the system. The dotted 
rectangle (1) refers to first flight test. In the language of block 
algebra, it represents the open loop of the airplane alone. The 
dotted rectangle (2) refers to the second flight test and represents 
the closed loop of the combination (aircraft plus pilot). The 
dotted rectangle (3) represents the open loop of the (aircraft plus 
pilot ) combination. 

The record of the first flight test (Fig. 1) provides data neces 
sary to compute the open loop transfer function of the airplane 
alone: 


F(jw) (1) 


The record of the second flight test (Fig. 2, starting from time 
ty) gives the data for computation of the transfer function: 
V(jw) (2) 
which refers to the closed loop of the (aircraft plus pilot) com 
bination. 
The data reduction that leads to the transfer function is a 
Where IBM machines are available, the 
Otherwise, the 


laborious process. 
method given in reference 4 is recommended. 
methods outlined in references 5 and 6 can be used. 

The block aigebra (Fig. 3) as given in the theory of servo- 
mechanisms” * gives the open loop of the (aircraft plus pilot) 
combination [see dotted rectangle (3)]: 

H(jw) = F(jw) X G(jw) = ao/(atrim — ao) (3) 
thus, 


G(jw) = ao/(atrim — ao) F(jw) (3a) 


Also, by the rules of the block algebra, the transfer function of the 
closed loop is 


V(jw) = H(jw)/(1 + H(jw)] (4) 
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From Eq. (4) 
H(jw) = Y(jw)/{1 — VGw)]) (5 


Eq. (5) substituted into Eq. (3) gives the transfer function of the 
human pilot as 


G(jw) = V(jw)/F(jw)[1 — V(jw)] (6 
Eq. (3a) or (6) can be used to compute the frequency response of 
the human pilot. By determining the real and imaginary part 
of G(jw), the amplitude and the angle of this transfer function can 
be obtained and tabulated for different values of the angular 
frequency w 

Next, the amplitude and angle of G(jw) can be plotted versus 
w (Fig. 4). Also, a polar plot can be used. 

No doubt a statistical approach is needed to arrive at some 
statistical conclusions or to get average data for a great number 
of pilots examined or for an average fighter or bomber type of ait- 
plane. 

The flight-test data can be supplemented using as a simulator 
a modern Link trainer, having provision for changing the damp- 
ing and sensitivity of the system. 

The above method gives, in the hand of the medical pro- 
fession, an exact tool to measure the behavior of the human pilot 
Cooperation with some engineering outfit in data reduction and 
program planning is, of course, desirable. 
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Compressible Laminar Boundary Layer on a 
Heat-Insulated Body* 


N. Rott 

Associate Professor of Aeronautical Engineering, Cornell University 
Ithaca, N.Y 

September 5, 1952 


ONSIDER A HEAT-INSULATED BODY submerged in a two-dimen 
C sional compressible gas flow. Suppose the isentropic flow 
outside the boundary layer is known—i.e., the velocity UU and 
Mach Number J\/ are given along the body contour. 


tropic temperature 7) outside the layer is 


To/T, = 1 + [(9 1)/2] M? (1) 


The isen 


(7) = stagnation temperature). The wall temperature 7°, over 
the surface portion where the boundary layer is laminar is given 
toa good approximation by the formula 

r 

T\ 


= 1+ 7(¢) |(y 1)/2] M? (2) 
where 7(o) is the “recovery factor,’’ which is a function of the 


Prandtl number o alone. For @ not too far from unity, 


r=Voe (0) 
and the connection between o and the adiabatic exponent y may 
be taken from Eucken’s approximate formula—viz., 

a0 = 4y/(9y — 5) (4) 
Theoretical support for Eq. (2) (originally derived by Pohl- 
hausen for the flat plate only) was given by calculations of Eckert 
and Drewitz! and Tifford,? carried out for velocity distributions 
U ~ xk. Experimental evidence was provided by Eckert and 
Weise* and by Ryan‘ through measurements on heat-insulated 
circular cylinders 
For the temperature 7 inside the boundary layer, at a point 
with local velocity “, an interpolation formula may be used which 


fulfills the boundary conditions 


7 ¥ l u? 
— = 1 + rao) M2 1 — : (5) 
7) 2 U? 


Eq. (5) is exact for o = r = 1, but for o ¥ 1 it gives a somewhat 


namely, 


less accurate result in general than for the special value at the 
wall. Here, it is proposed to accept Eq. (5) as an approximate 
solution for the “thermal”? part of the boundary-layer problem 
with the purpose of computing the density p and viscosity wu by 
part of the prob- 


Eq. (5), so that subsequently the “momentum 
lem may be attacked. 
The density variation across the boundary layer is given by the 
law 
pI = pl, = Sal « (6) 
As to the variation of the viscosity, a further simplification will 
be needed, introduced first by Chapman and Rubesin.® Let us 
assume a variation of uw proportional to 7 across the boundary 


*This research was supported by the U.S.A.F. under Contract No 
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layer, 

b/w = T/Ty (7) 
but let the wall value nu» (7) be computed by any more accurate 
formula—e.g., 


w/o = (T~/To)® (8) 


with 0.5 < w < 1.5—or by Sutherland’s law. Naturally, some 
inconsistency is introduced by evaluating « according to different 
laws along and across the boundary layer, but it is felt that it is 
more important to calculate 4» accurately, whereas the variation 
across the layer is less significant, especially for the heat-insulated 
body with vanishing normal temperature gradient on the surface. 

The purpose of the present note is to point out that, once the 
approximations stated above are accepted, the remaining prob- 
lem may be reduced exactly to an incompressible laminar bound- 
ary-layer problem by a slight extension of a transformation given 
previously by Illingworth® and Stewartson.’ The extended trans- 
formation will be given without proof, since its verification closely 
follows the derivation of the original transformation. The prob- 
lem is reduced to the incompressible case by introduction of new 
coordinates x;, 2:;, parallel and perpendicular to the surface, which 


are connected with the actual coordinates x, s by the equations 


wx - ly (y 1)]) + (r/2) 7 
_ (#) | oe (9) 
0 To T wuo 


T,\0/(7 — DI + (7/2 “s 
5 = ( ) = ds (10) 
T; Oo Ppl 


and in the incompressible case, a velocity LU’; has to be taken out- 
side the layer which is connected with LU’ by the formula 
U; = (T/T)"/? U (11) 
The original Illingworth-Stewartson transformation, which was 
derived for ¢o = r = w = 1, is included; in this special case, the 
equations for the “thermal” part of the solution are exact, and the 
whole solution is also exact. The present extension allows one to 
eliminate the restrictions concerning the material properties of the 
gas, naturally at the price of the approximations discussed above. 
It is possible to apply the transformation, Eqs. (9) to (11), to 
the “quadrature”? formulas of incompressible boundary-layer 
theory, which are simplified versions of the Pohlhausen method 
This was carried through for the original Illingworth-Stewartson 
transformation in a recent paper by Rott and Crabtree,’ so that, 
here again, only the results will be presented 
The momentum thickness @ of the compressible boundary 


layer is given by the formula 


i [2/(y 1)] 2) 
6? = 0.45 » (2) U-* X 
Ty 


r,t - BU - & Ty, 
(') HY Usdx (12) 
0 To TP wpo 


or, for y = 1.4 and o = 0.737 from Eq. (4), and from Eq. (8) 


T, \3-32 
e2 = 0) ton (7) Us x 
T; 
“* (7,\1-82/ 7, \l - & 
(7) (Z U'dx (12a) 
0 T 


Here, vp is the kinematic viscosity at the stagnation point 
The (modified) Pohlhausen parameter m; (see reference 8) for 
the corresponding incompressible flow is 
dU @ ( T,\8 — 20/7 — D Tyg 
== om ies (13) 
dx w\T»o To20 


which becomes, with the numerical values as before, 


ave (ny 7... 9 
n= 
dx w\T% To 


(13a) 
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The shape parameter /7; and the non-dimensional shearing stress 
J; are known functions of m; (see reference 8). In the compres- 


sible flow, H follows from the relation 
H = (T./1T;)Hi + (Tw/TM) — 1 (14) 
and the shearing stress is given by the formula 
tw = li ww(T1/Tw) (U/0) (15) 


As a simple check, consider the flat plate, VU = constant. For 
the coefficient of friction, c;, it is easily found from Eqs. (12), 
(15), and (8) that 

cy = cy, (T1/T)| — ”/2 (16) 
where cy, is the skin friction coefficient calculated with values 
p, and yw; outside the boundary layer by the formula for incom- 
pressible flow. Eq. (16) contains the exact result that, forw = 1 
and arbitrary o, cy = Cy,. 

If w # 1 but o = 1, comparison with exact calculations of 
Hantzsche and Wendt shows that Eq. (16) still reproduces the 
main part of the compressibility effect on c;, although it is not 
exact. However, even for moderately supersonic Mach Numbers, 
the difference is completely negligible. 
co, Eq. (16) gives cy-values that are 6 and 15 per cent too low 


In the limiting case J —> 


In the general case, when 


w * land o ¥ 1, comparison with Crocco’s calculations!® 


for w = 0.8 and w = 0.5, respectively. 
shows 
differences of similar character. 

It is believed that for not too high Mach Numbers and for a 
heat-insulated body (which is an idealized case anyway ) the pres- 


Potential Limit in Transonic Flow 


Pei Chi Chou 
Aerodynamicist, Republic Aviation Corporation, New York 
September 4, 1952 


N REFERENCE I, A CRITERION OF COLLAPSE of potential flow is 

formulated for flow with a local supersonic region bounded 
by the airfoil and the line of transition from subsonic to super- 
sonic velocity. This criterion is in terms of the local velocity and 
local airfoil curvature. Nothing is mentioned about the free- 
stream velocity and airfoil 
note, following the method given in reference 2, numerical 
solutions are found for flow past different families of airfoils. The 
relations among the limiting free-stream Mach Number, limiting 
maximum local Mach Number, the airfoil thickness ratio, and the 
local airfoil curvature are plotted explicitly. The result shows 
agreement with reference 1. 

In reference 2, an example has been worked out for a Kaplan 


overall configuration. In this 


section by a simplified iteration process. For comparison, two 
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ent simple approach is sufficiently accurate. Extension for bod. 


ies of revolution is easily obtained, as was shown in reference 8. 
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modified Kaplan sections are solved by the same method 
here. The parametric representation of these modified sections 
are: 

Modified Section I 


x = cost 


| 


y = (7/5) (4 sin ¢ — sin 3¢) 
Modified Section II 


x = cost 


y = (7/6) (5 sin t — sin 3¢) 


The configuration of these sections is shown in Fig. 1. The 
curvatures C at midchord are 37, 2.67, and 2.337 for Kaplan sec- 
tion, the modified section I, and the modified section II, respec- 
tively, where 7 is the thickness ratio. The general forms of the 
three are alike, except that the last two differ from the first in that 
C is not zero at leading edge and, hence, introduces stagnation 
points. The simplified equation used here, however, is valid only 
within a rectangle of a length of about 0.3 chord at the midchord 


Therefore the effect of stagnation point can be neglected. 













— KAPLAN 


——— MODIFIED I 
a MODIFIED 1 


FIG. 1 K APLAN 










AND MODIFIED KAPLAN SECTIONS 








ING MAX. LOCAL MACH NO. 
= a 3 Ss 
r Pad 


Min, LIMITING 
S 
~~ 





The r 
approxi 


Modi} 


(Or) m 
Modi 


( Pr) n 


Let t 
to the | 
waves 
tion I 
limiting 
local M 
vs. Moz 
equal tk 
conclus 

(1) F 
ness rat 
and dec 
can sus 
stream 
sidering 
referenc 

(2) If 
foil mic 
some d 

(3) I 
ratio if 

The 
than tk 
ployed. 
descrip 


1 Nike 
Region a 
May, 19 





»d- 


b der 
hrt- 











READERS’ 














. \—-— KAPLAN 
S |—-— MODIFIED I 
: —~— MODIFIED I 
|----- -CONST. C 
a | 
< mea 5 
Y |.08 ; <—e 
5 
a | 
s! | 
v ae > 
z . > 
~@ | ZN 
= 
5 102 ee 
F | 
| ee ee 

















0.8 09 


0.7 
Mo. LIMITING FREE STREAM MACH NO. 


FIG. 2 Mp, VS Mog, CURVES 


The resulting maximum velocity increments for the first three 
approximations are: 
Modified Section I 


(¢r)maz. = (7/8) (1.4000 + 0.7153K + 1.0273K? +...) 


Modified Section II 


(dz)mar. = (7/B) (1.3333 + 0.6517K + 0.9220K? + ...) 
Let the subscript ‘““’’ denote the limiting value corresponding 

to the breakdown of potential flow or the occurrence of shock 

The estimated Ky, are 0.375 and 0.38 for modified sec 


The corresponding values of 


waves. 
tion I and modified section II. 
limiting free-stream Mach Number MJ, and limiting maximum 
local Mach Number Mj, are plotted in Fig. 2. Adding the Mnz 
vs. Moz curve of Kaplan section to Fig. 2 and connecting the 
equal thickness ratio lines and equal curvature lines, the following 
conclusions are drawn: 

(1) For the same family of airfoil sections, an increase in thick- 
ness ratio (hence, local curvature) results in an increase of Mz 
and decrease of M/o,. In other words, a thicker airfoil in a family 
can sustain a higher local Mach Number but at a lower free 
stream Mach number without introducing shock waves. Con- 
sidering M1 and local airfoil curvature alone, this agrees with 
reference 1. 

(2) If the thickness ratio is kept constant, an increase in the air- 
foil midchord curvature raises the M,,, considerably at a cost of 
some decrease of Moz 

(3) Both Moz and Mz increase with a decrease in thickness 
ratio if the midchord curvature can be kept unchanged. 

The values of M,,, from these curves are considerably lower 
than the actual values, since only three approximations are em 
ployed. However, the above conclusions give a good qualitative 
description. 
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Note on the Dynamic Stability of a Missile in 
Rolling Flight 


Arthur E. Bryson, Jr. 

Guided Missile Laboratories, Hughes Aircraft Company, 
Culver City, Calif. 

September 26, 1952 


T° Ray E. Bouz’s PAPER “Dynamic Stability of a Missile in 
Rolling Flight’”’ (JouRNAL OF THE AERONAUTICAL SCIENCES, 
June, 1952), we believe there are several errors that should be 
corrected. They are enumerated below in the order in which 
they appeared in the subject paper. 

(1) In Egs. (8e) and (3f) the quantity called P; can only be the 
instantaneous roll rate of the missile. Eq. (4) is confusing, since 
it implies that P; is made up of the instantaneous roll rate plus 
“an angular velocity w, about the XY axis due to precessional mo- 
tion.” 

(2) The inclusion of forces and moments due to misaligned and 
canted surfaces and jet misalignment in a linearized stability 
analysis is obviously unnecessary, since their effect is simply to 
change the trim attitude of the missile. 

(3) If one considers the small forces due to pitching and yawing 
angular velocities, as Mr. Bolz does, to be consistent one should 
also include the Magnus forces and moments due to pitching and 
yawing angular velocities. Since none of these quantities affects 
the stability to any important degree, they can be neglected. 

(4) In Table II, item 4, it is stated that Kn, is always positive 
We believe it is almost always negative for the manner in which 
it is defined in the paper. 

(5) In Table I, item 7, Kp should be equal to Cp (.S/2d?). Also, 
22), should read 


the term RK. — 2Kp, which occurs in Eq. (22), 


RK. — Kp. This is the lift coefficient curve slope, which is always 
positive neat a = 0. Likewise, the term 2K in Eq. (22) should 
read K. 


(6) The third stability criterion (III) on p. 400 is already im 
plied by the other two stability conditions and, hence, is unneces- 
sary. 

(7) The transformation from time to distance along the flight 
path in Eqs. (15) is just as valid for subsonic as supersonic veloci 
ties. Thus the conclusion under (H) on p. 401—that dynamic 
stability is independent of velocity (except for variation of the 
coefficients with Mach Number)—applies equally well to sub- 
sonic velocities as it does to supersonic velocities 

(8) If reference axes are used which pitch and yaw with the mis 
sile but do not roll with it, the stability analysis is simplified and 
clarified (e.g., see reference 3 of Mr. Bolz’s paper). If ¢ and 7 
are the complex angular and linear velocities instantaneously re- 


solved on this nonrolling reference system, then 


where ¢g and 7 are Mr. Bolz’s quantities. Taking the derivative 
with respect to time, 
¢=e($+iPe); 7 


Thus several terms could be eliminated from Eqs. (8a) and (8b) of 


= Fania (1 + iPn) 


Mr. Bolz’s paper. 

(9) In Tables I and II, item 3, reference is made to forces and 
moments ‘due to rolling combined with pitching and yawing 
Apparently it was intended to read “. . . com- 
These forces and mo- 


angular velocities.’ 
bined with angles of attack and yaw.” 
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ments are actually the same as those under item 4; in nonrolling 
coordinates, one would write the complex moment due to accel- 
erated motion as 


Kn, 03% 
Making the change to the rolling coordinates of Bolz, 7 is replaced 
by 
e **(d/dt) (ne?*) = 4 + iPn 


and, hence, the complex moment due to accelerated motion in 
Bolz’s rolling coordinates is given by 

Km gpd + iPr) 
Therefore, Kz, = —K,. and Kmp = —K,,.. in Bolz’s notation. 
Since Kz, has been neglected, Kz, should also be neglected. 

(10) We find the ‘‘discussion of Stability” rather unenlightening 
and in some cases misleading. For example, a missile that is 
dynamically stable (and, hence, also statically stable) for zero roll 
rate cannot be made dynamically unstable by the presence of a 
roll rate except under rather unusual conditions. This is shown in 
the following analysis, which is similar to that of Bolz except that 
we believe all errors have been eliminated, unessential forces and 
moments have been neglected, a simpler reference system has 
been used, and dimensionless quantities have been used through- 


out. 
SYMBOLS 
a = —i[(v + tw)/V], complex angle of attack 
w = q + ir, complex pitching velocity 
l = 2m/pS, ballistic length 
YA = VA m, radius of gyration in roll 
rapa OU /B/m, radius of gyration in pitch (yaw) 
Ce = Ci, — Cp, lift coefficient curve slope 
Cine = Magnus lift due to angle of attack + [(pV?/2)S(Pd/- 
V)a] 
Cu, = moment coefficient curve slope 
Cu, = 0Cy/0(ad/V), damping in plunge 
CM pa = Magnus moment due to angle of attack + [(pV?/2) x 


Sd(Pd/V)a] 


Cu, = 0Cy/O(qd/V), damping in pitch 


and all other symbols not defined here are the same as Bolz’s 
symbols or are defined in the text. All lift coefficient curve 
slopes are always positive; the moment coefficient curve slopes 
are positive if the point of action of the force creating the moment 
is ahead of the center of gravity. 
ANALYSIS 

The two complex equations of motion of a coasting missile (zero 

thrust) referred to axes that pitch and yaw with the missile but 


do not roll with it are 


— 2 Pd 
o-a- i Ci, — iy, Vy a=Q0 


: A P V*d C wd £C ad ys (1) 
3 a w ire? My y Mey 


. sel Pd 
Cu, — 1¢ M nw V al=0 
' V(r) 
z= [ se dr 
A d 


distance along the flight path in diameters, ( )’ = derivative with 


respect to z, and 
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Ci = (d/l)CL 

Cu = (d/l) (d/rg)? Cu 
CMa = d/Kd/ra)?Cat yg 
Py) = (A/B)(Pd/V) 


Then Eqs. (1) become 


d 
V 


w-—a’— (Ci, _ ICL ygPoa = 0 
* a ee 
V [w iin (Cm, a 1Po)w] = 


; = = 
Cm, a (Cu, iC Poe = () 


It will be noticed that all imaginary terms (the Magnus force 
and moment and the gyroscopic torque) are proportional to the 
roll rate. 

The characteristic equation of Eqs. (2), divided by d/ V, is 


d+ Cr, — iCr, Po 


i] 


A- (Cu, + iPo) Cur, -+- Cu, - Cu, 


or 
A? + (a; + iPoaz)A + bi + 1Pobe = O (3 
where 


ay = x = (Cu, = Cu) 


a2 = —1+ Cpa 
by = a — Cu, Cia — PPC ie 
be = Cung + Cigg — Chg 


The conjugate of Eq. (3) is, of course, also a characteristic 
equation. The product of Eq. (3) with its conjugate yields the 
stability quartic with real coefficients, which would have been ob 
tained by the usual analysis (i.e., not using complex notation). 

For stability, the real part of the roots of Eq. (3) must be neg- 
ative. Writing down this condition in a manner similar to Bolz, 
the two stability conditions become 


a, > 0 
by>be/ay (bo/a; = @ )Po? 


Substituting in the expressions for a, d2, b;, and b., the stability 


conditions become 
= (Cug + Cu ;) > 0 


—Cyw. > 


a 
(cs. = Choa = CM pq) (Cmq + Cu, — Chow ~ Cas - 
ieee: = o 
(Cx, = (Car, + Cm,)]? 
where products and squares of the barred coefficients have been 
neglected with respect to the barred coefficients themselves 
(which, in magnitude, are all very much less than unity). 
If we let 
1 = —(Cu, + Cm ,) CL, 


a = —(Crpa + CM nq) CL, 


then, finally, the stability conditions read 


5s, > -—1 (4a) 
(1 + Se) (s1 — Se) 
BM, < ote = carp (4b) 
: (1 + s:)? 
where M, = Cu(p V? 2)Sd = pitching moment curve slope. 


Eq. (4a) states that there must be positive aerodynamic damp- 
ing. Eq. (4b) is a condition on the static stability; if P = 0, it 
says that the pitching moment curve slope must be negative. 
Now the only possibility that a statically stable missile (Ma < 
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0) with positive aerodynamic damping (s, > —1) could be Incidentally, if we let s = [(L + se)/(1 + s:)] — (1/2), we 
made dynamically unstable by the presence of a roll rate is for may rewrite Eq. (4b) as 

1 $2) (Ss) — Se) < O 

(1 + se) (1 — Se BM,/(AP)? < (1/4) — s? (5) 


For this to take place either sp < —1 or sz > s;. Under these 
conditions a statically unstable missile (M/q > 0) with positive 
aerodynamic damping (s; > —1) could not be spin-stabilized by 
Upon examining these conditions it 


showing that the right-hand side has a stationary value in the 
neighborhood of s = 0 whichisa maximum. Eq. (5) withs = Ois 
precisely the condition for stability of a heavy symmetrical top 
with one point fixed; it is also the approximate stability condition 
that is often used for spin-stabilized projectiles, and it can be 
seen from Eq. (5) that this is a very good approximation, since 
usually |s|< 1. 


increasing the roll rate. 
would appear to us that the missile would have to be of a peculiar 
design for them to be satisfied, since the Magnus force and mo 
ment are usually so small that 52 <_1 and for almost all mis- 








siles and projectiles s; ~ 1. 
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